analysis_INTT_CW/DST_MC/INTTEta.h

0001 #ifndef INTTEta_h
0002 #define INTTEta_h
0003
0004 #include "INTTXYvtxEvt.h"
0005 #include "MegaTrackFinder.h"
0006 #include "ana_map_folder/ana_map_v1.h"
0007
0008 // todo : if the centrality bin or the eta and z regions are changed, modify the "ana_map_v1.h" first,
0009 // todo : and change the following "namespace" name
0010 namespace ana_map_version = ANA_MAP_V3;
0011
0012 class INTTEta : public INTTXYvtxEvt
0013 {
0014     public:
0015         INTTEta(
0016             string run_type,
0017             string out_folder_directory,
0018             pair<double,double> beam_origin,
0019             int geo_mode_id,
0020             double phi_diff_cut = 0.11,
0021             pair<double, double> DCA_cut = {-1,1},
0022             int N_clu_cutl = 20,
0023             int N_clu_cut = 10000,
0024             bool draw_event_display = true,
0025             double peek = 3.32405,
0026             double angle_diff_new_l = 0.0,
0027             double angle_diff_new_r = 3,
0028             bool print_message_opt = true
0029         ):
0030         INTTXYvtxEvt(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, 2.5, 9)
0031         {
0032             track_cluster_ratio_1D.clear();
0033             track_cluster_ratio_1D_MC.clear();
0034             dNdeta_1D.clear();
0035             track_eta_phi_2D.clear();
0036             track_eta_z_2D.clear();
0037             dNdeta_1D_MC.clear();
0038             dNdeta_1D_MC_edge_eta_cut.clear();
0039             final_dNdeta_1D.clear();
0040             final_dNdeta_1D_MC.clear();
0041             final_dNdeta_multi_1D.clear();
0042             track_eta_z_2D_MC.clear();
0043             track_delta_phi_1D.clear();
0044             track_DCA_distance.clear();
0045             final_track_delta_phi_1D.clear();
0046             final_track_multi_delta_phi_1D.clear();
0047
0048             track_delta_eta_1D.clear();
0049             track_delta_eta_1D_post.clear();
0050             track_phi_DCA_2D.clear();
0051             track_DeltaPhi_eta_2D.clear();
0052             track_DeltaPhi_DeltaEta_2D.clear();
0053             track_ratio_1D.clear();
0054
0055             proto_pair_index.clear();
0056             proto_pair_delta_phi_abs.clear();
0057             proto_pair_delta_phi.clear();
0058             inner_used_clu.clear();
0059             outer_used_clu.clear();
0060
0061             coarse_eta_z_2D_MC.clear();
0062             coarse_eta_z_2D_fulleta_MC.clear();
0063             coarse_Reco_SignalNTracklet_Single_eta_z_2D.clear();
0064             coarse_Reco_SignalNTracklet_Multi_eta_z_2D.clear();
0065
0066             out_recotrack_eta_d.clear();
0067             out_recotrack_phi_d.clear();
0068             out_truetrack_eta_d.clear();
0069             out_truetrack_phi_d.clear();
0070
0071             out_track_eta_i.clear();
0072             out_track_delta_phi_d.clear();
0073             clu_multi_used_tight = 0;
0074             clu_multi_used_loose = 0;
0075             effective_total_NClus = 0;
0076
0077             out_N2Clu_track = 0;
0078             out_N3Clu_track = 0;
0079             out_N4Clu_track = 0;
0080
0081             eta_z_convert_map.clear();
0082             eta_z_convert_inverse_map.clear();
0083             eta_z_convert_map_index = 0;
0084
0085             // std::fill(std::begin(inner_clu_phi_map), std::end(inner_clu_phi_map), std::vector<pair<bool,clu_info>>());
0086             // std::fill(std::begin(outer_clu_phi_map), std::end(outer_clu_phi_map), std::vector<pair<bool,clu_info>>());
0087             inner_clu_phi_map.clear();
0088             outer_clu_phi_map.clear();
0089             inner_clu_phi_map = vector<vector<pair<bool,clu_info>>>(360);
0090             outer_clu_phi_map = vector<vector<pair<bool,clu_info>>>(360);
0091
0092             mega_track_finder = new MegaTrackFinder(run_type, out_folder_directory, centrality_map.size(), beam_origin, true); // note : "true" means we marked the clusters as used
0093
0094             // note : 10 is for the centrality bin
0095             N_GoodEvent_vec = vector<long long>(10,0); // todo : if the centrality bin is changed, the vector should be updated
0096
0097             InitHist();
0098             InitGraph();
0099             InitTree();
0100
0101             cout<<"test, eta_region.size() : "<<eta_region.size()<<endl;
0102             cout<<"test, z_region.size() : "<<z_region.size()<<endl;
0103             for (int i = 0; i < z_region.size() - 1; i++) // note : Y axis in the TH2F
0104             {
0105                 for (int j = 0; j < eta_region.size() - 1; j++) // note : X axis in the TH2F
0106                 {
0107                     eta_z_convert_map[std::to_string(j+1) + "_" + std::to_string(i+1)] = eta_z_convert_map_index;
0108                     eta_z_convert_inverse_map[eta_z_convert_map_index] = {j+1, i+1};
0109                     eta_z_convert_map_index += 1;
0110                 }
0111             }
0112
0113             cout<<"In INTTEta, centrality_region_size: "<<centrality_region_size<<" eta_region_size: "<<eta_region_size<<" z_region_size: "<<z_region.size() - 1<<endl;
0114             good_tracklet_counting = std::vector<std::vector<int>>(centrality_region_size, std::vector<int>(((eta_region.size() - 1) * (z_region.size() - 1)), 0));
0115             good_tracklet_multi_counting = std::vector<std::vector<int>>(centrality_region_size, std::vector<int>(((eta_region.size() - 1) * (z_region.size() - 1)), 0));
0116
0117             out_folder_directory_evt = out_folder_directory + "/evt_display";
0118             system(Form("mkdir %s", out_folder_directory_evt.c_str()));
0119
0120             loop_NGoodPair = 0;
0121             evt_NTrack = 0;
0122             evt_NTrack_MC = 0;
0123             N_GoodEvent = 0;
0124             return;
0125         };
0126
0127         INTTEta(string run_type, string out_folder_directory):INTTXYvtxEvt(run_type, out_folder_directory)
0128         {
0129             return;
0130         };
0131
0132         virtual void ProcessEvt(
0133             int event_i,
0134             vector<clu_info> temp_sPH_inner_nocolumn_vec,
0135             vector<clu_info> temp_sPH_outer_nocolumn_vec,
0136             vector<vector<double>> temp_sPH_nocolumn_vec,
0137             vector<vector<double>> temp_sPH_nocolumn_rz_vec,
0138             int NvtxMC,
0139             vector<double> TrigvtxMC,
0140             bool PhiCheckTag,
0141             Long64_t bco_full,
0142             pair<double,double> evt_z,
0143             int centrality_bin,
0144             vector<vector<float>> true_track_info,
0145             double zvtx_weighting
0146         );
0147         // void ProcessEvtMC(int event_i, vector<float> true_track_info, vector<double> TrigvtxMC);
0148         void ClearEvt() override;
0149         void PrintPlots() override;
0150         void EndRun() override;
0151
0152     protected:
0153         // note : {5, 15, 25, 35, 45, 55, 65, 75, 85, 95, inclusive}
0154         vector<TH1F *> track_cluster_ratio_1D;
0155         vector<TH1F *> track_cluster_ratio_1D_MC;
0156         vector<TH1F *> dNdeta_1D;
0157         vector<TH2F *> track_eta_phi_2D;
0158         vector<TH2F *> track_eta_z_2D;
0159
0160         vector<TH1F *> dNdeta_1D_MC;
0161         vector<TH1F *> dNdeta_1D_MC_edge_eta_cut;
0162         vector<TH2F *> track_eta_z_2D_MC;
0163
0164         vector<TH1F *> track_delta_eta_1D;
0165         vector<TH1F *> track_delta_eta_1D_post;
0166         vector<TH1F *> track_delta_phi_1D;
0167         vector<TH1F *> track_DCA_distance;
0168         vector<TH2F *> track_phi_DCA_2D;
0169
0170         vector<TH2F *> track_DeltaPhi_eta_2D;
0171         vector<TH2F *> track_DeltaPhi_DeltaEta_2D;
0172
0173         TH2F * track_correlation_2D;
0174         TH2F * track_ratio_2D;
0175         vector<TH1F *> track_ratio_1D; // note : N reco track /N true track
0176
0177         TH2F * reco_eta_correlation_2D;
0178         TH2F * reco_eta_diff_reco3P_2D;
0179         TH1F * reco_eta_diff_1D;
0180
0181         TH2F * clu_used_centrality_2D;
0182
0183         vector<vector<TH1F *>> final_track_delta_phi_1D;       // note : centrality_bin and different eta
0184         vector<vector<TH1F *>> final_track_multi_delta_phi_1D; // note : centrality_bin and different eta
0185         // vector<vector<int>> good_tracklet_counting;
0186         vector<TH1F *> final_dNdeta_1D;
0187         vector<TH1F *> final_dNdeta_multi_1D;
0188         vector<TH1F *> final_dNdeta_1D_MC;
0189
0190         TGraph * evt_reco_track_gr_All; // note : all the tracklets in one event, no selection
0191         TGraph * evt_reco_track_gr_PhiLoose; // note : the tracklets that pass the loose phi selection
0192         TGraph * evt_reco_track_gr_Z; // note : the tracklets that pass the z selection and the loose phi selection
0193         TGraph * evt_reco_track_gr_ZDCA; // note : the tracklets that pass the z and phi selection, and the loose phi selection
0194         TGraph * evt_reco_track_gr_ZDCAPhi; // note : the tracklets that pass the z, phi and DCA selection, and the loose phi selection
0195         TGraph * evt_true_track_gr; // note : the true tracklets in one event, no selection
0196
0197         MegaTrackFinder * mega_track_finder; // note : the class that handles the 3/4-cluster tracklet finder
0198         TH2F * NClu3_track_centrality_2D;
0199         TH2F * NClu4_track_centrality_2D;
0200         TH1F * cluster4_track_phi_1D;
0201         TH1F * cluster3_track_phi_1D;
0202         TH1F * cluster3_inner_track_eta_1D;
0203         TH1F * cluster3_inner_track_phi_1D;
0204         TH1F * cluster3_outer_track_eta_1D;
0205         TH1F * cluster3_outer_track_phi_1D;
0206         TH1F * clu4_track_ReducedChi2_1D;
0207         TH1F * clu3_track_ReducedChi2_1D;
0208         TH1F * mega_track_eta_1D;
0209         TH2F * mega_track_finding_ratio_2D;
0210         TH2F * mega_track_contamination_2D;
0211
0212         // note : diagnostic plots
0213         TH1F * check_inner_layer_clu_phi_1D;
0214         TH1F * check_outer_layer_clu_phi_1D;
0215
0216         // note : for the comparison with other files, MC, for example
0217         // note : after all the selections, so it is exclusive
0218         // note : cut : cluster size, cluster adc
0219         // note : min bias cut
0220         // note : zvtx quality cut
0221         // note: zvtx acceptance cut
0222         TH1F * exclusive_NClus_inner;
0223         TH1F * exclusive_NClus_outer;
0224         TH1F * exclusive_NClus_sum;
0225
0226         TH1F * exclusive_cluster_inner_adc;
0227         TH1F * exclusive_cluster_outer_adc;
0228
0229         TH2F * exclusive_cluster_inner_eta_adc_2D;
0230         TH2F * exclusive_cluster_outer_eta_adc_2D;
0231
0232         TH2F * exclusive_cluster_inner_eta_phi_2D;
0233         TH2F * exclusive_cluster_outer_eta_phi_2D;
0234
0235         TH1F * exclusive_cluster_inner_eta;
0236         TH1F * exclusive_cluster_inner_phi;
0237         TH1F * exclusive_cluster_outer_eta;
0238         TH1F * exclusive_cluster_outer_phi;
0239         TH1F * exclusive_cluster_all_eta;
0240         TH1F * exclusive_cluster_all_phi;
0241
0242         TH1F * exclusive_tight_tracklet_eta; // note : signal region only
0243         TH1F * exclusive_tight_tracklet_phi; // note : signal region only
0244         TH1F * exclusive_loose_tracklet_eta; // note : signal region only
0245         TH1F * exclusive_loose_tracklet_phi; // note : signal region only
0246
0247
0248         // vector<pair<bool,clu_info>> inner_clu_phi_map[360]; // note: phi
0249         // vector<pair<bool,clu_info>> outer_clu_phi_map[360]; // note: phi
0250         vector<vector<pair<bool,clu_info>>> inner_clu_phi_map; // note: phi
0251         vector<vector<pair<bool,clu_info>>> outer_clu_phi_map; // note: phi
0252
0253         string out_folder_directory_evt;
0254
0255         int draw_evt_cut = 15; // todo : change the cut value (number of true track)
0256         int print_plot_ratio = 3; // note : print the event display every 3 events
0257         int inner_NClu;
0258         int outer_NClu;
0259         int loop_NGoodPair;
0260         int evt_NTrack;
0261         int evt_NTrack_MC;
0262         long long N_GoodEvent;
0263         int clu_multi_used_loose;
0264         int clu_multi_used_tight;
0265         int effective_total_NClus;
0266         double pair_delta_phi;
0267         pair<int,int> pair_outer_index;
0268         // double INTT_layer_R = 71.88; // note: the innermost, B0L0, layer radius
0269         double INTT_layer_R[4] = {71.88, 77.32, 96.8, 102.62}; // note : the radii of the INTT layers
0270
0271         vector<long long> N_GoodEvent_vec;
0272         pair<double,double> Get_eta_pair;
0273         vector<vector<double>> final_eta_entry; // note : not used
0274
0275
0276         vector<vector<int>> proto_pair_index;
0277         vector<vector<double>> proto_pair_delta_phi;
0278         vector<double> proto_pair_delta_phi_abs;
0279         map<string,int> inner_used_clu;
0280         map<string,int> outer_used_clu;
0281
0282         double Clus_InnerPhi_Offset_1;
0283         double Clus_InnerPhi_Offset_2;
0284         double Clus_OuterPhi_Offset_1;
0285         double Clus_OuterPhi_Offset_2;
0286
0287         double N_reco_cluster_short = 0;
0288
0289
0290         TH2F * track_cluster_ratio_multiplicity_2D; // note : x : NClus / NTracks, y : ratio
0291         TH2F * track_cluster_ratio_multiplicity_2D_MC; // note : x : NClus / NTracks, y : ratio
0292         TGraphErrors * track_gr;
0293
0294         TFile * out_file;
0295         TTree * tree_out;
0296         int out_eID;
0297         int out_evt_centrality_bin;
0298         int out_NTrueTrack;
0299         int out_total_NClus;
0300         int out_N2Clu_track; // note : method tight
0301         int out_N3Clu_track; // note : mega cluster track 3-cluster track
0302         int out_N4Clu_track; // note : mega cluster track 4-cluster track
0303         double out_evt_zvtx;
0304         double out_true_zvtx;
0305         vector<double> out_recotrack_eta_d;
0306         vector<double> out_recotrack_phi_d;
0307         vector<double> out_truetrack_eta_d;
0308         vector<double> out_truetrack_phi_d;
0309         vector<int> out_track_eta_i;
0310         vector<double> out_track_delta_phi_d;
0311
0312         TH1F * eta_region_hist;
0313
0314         map<int,int> centrality_map = ana_map_version::centrality_map;
0315         vector<string> centrality_region = ana_map_version::centrality_region;
0316         // note : the following two vectors tell the region of the eta and z, including the edges in both sides.
0317         // note : so when talking about the number of bins, we have to do size() - 1
0318         vector<double> eta_region = ana_map_version::eta_region;
0319         vector<double> z_region = ana_map_version::z_region;
0320         pair<double, double> selected_z_region_id = ana_map_version::selected_z_region_id;
0321         double signal_region = ana_map_version::signal_region; // note : the signal region of the delta phi distribution, unit : degree
0322
0323         // note : when filling in a z-eta pair, TH2F returns a int, we have a function to convert that number into the index position in X and Y
0324         // note : Once we have the position index in X and Y, we then convert the X and Y into the index for 1D array
0325         int eta_z_convert_map_index;
0326         map<string,int> eta_z_convert_map;
0327         map<int, pair<int,int>> eta_z_convert_inverse_map;
0328
0329         TH2F * coarse_eta_z_map; // note : the reference of the eta and z binning setting
0330         vector<TH2F *> coarse_eta_z_2D_MC; // note : keep the number of true tracks in each eta and z bin, for each centrality
0331         vector<TH2F *> coarse_eta_z_2D_fulleta_MC; // note : no selection on the true level tracks
0332         vector<TH2F *> coarse_Reco_SignalNTracklet_Single_eta_z_2D; // note : keep the counting of the N reco tracklet in the signal region for each eta, z and centrality bin. (for the single-cluster-used tracklet)
0333         vector<TH2F *> coarse_Reco_SignalNTracklet_Multi_eta_z_2D;  // note : // note : keep the counting of the N reco tracklet in the signal region for each eta, z and centrality bin. (for the multi-cluster-used tracklet)
0334         TH2F * MBin_Z_evt_hist_2D;    // note : keep the number of event as a function of the z vertex and centrality bin
0335         TH2F * MBin_Z_evt_hist_2D_MC; // note : keep the number of event as a function of the z vertex and centrality bin, but the zvtx info. is given by the MC, the true zvtx
0336
0337
0338
0339         // todo : change the tight zvtx cut here
0340         int tight_zvtx_bin = 6; // todo : it can be run by run dependent // note : select the bin that you want to have the quick result out
0341
0342         const int centrality_region_size = centrality_region.size();
0343         const int  eta_region_size = eta_region.size() - 1;
0344         std::vector<std::vector<int>> good_tracklet_counting;
0345         std::vector<std::vector<int>> good_tracklet_multi_counting;
0346
0347         void InitHist();
0348         void InitGraph();
0349         void InitTree();
0350
0351         double grEY_stddev(TGraphErrors * input_grr);
0352         pair<double, double> mirrorPolynomial(double a, double b);
0353         pair<double,double> Get_eta(vector<double>p0, vector<double>p1, vector<double>p2);
0354         double convertTo360(double radian);
0355         void print_evt_plot(int event_i, int NTrueTrack, int innerNClu, int outerNClu);
0356         double get_clu_eta(vector<double> vertex, vector<double> clu_pos);
0357         // double get_dist_offset(TH1F * hist_in, int check_N_bin);
0358         double get_delta_phi(double angle_1, double angle_2);
0359         double get_track_phi(double inner_clu_phi_in, double delta_phi_in);
0360         pair<int,int> GetTH2BinXY(int histNbinsX, int histNbinsY, int binglobal);
0361         double GetTH2Index1D(pair<int,int> XY_index, int histNbinsX);
0362         void DrawEtaZGrid();
0363         double get_TH1F_Entries(TH1F * hist_in);
0364 };
0365
0366 void INTTEta::InitHist()
0367 {
0368     // double Eta_NBin = 145;
0369     // double Eta_Min = -2.9;
0370     // double Eta_Max = 2.9;
0371
0372     // todo: change the fine eta binning here
0373     double Eta_NBin = 61;
0374     double Eta_Min = -3.05;
0375     double Eta_Max = 3.05;
0376
0377     final_track_delta_phi_1D = vector<vector<TH1F *>>(centrality_region.size());
0378     final_track_multi_delta_phi_1D = vector<vector<TH1F *>>(centrality_region.size());
0379     final_eta_entry = vector<vector<double>>(centrality_region.size());
0380
0381     coarse_eta_z_map = new TH2F("","", eta_region.size() - 1, &eta_region[0], z_region.size() - 1, &z_region[0]);
0382     coarse_eta_z_map -> GetXaxis() -> SetTitle("#eta");
0383     coarse_eta_z_map -> GetYaxis() -> SetTitle("Z vertex [mm]");
0384     coarse_eta_z_map -> GetXaxis() -> SetNdivisions(505);
0385
0386     for (int i = 0; i < centrality_region.size(); i++)
0387     {
0388
0389         coarse_eta_z_2D_MC.push_back( new TH2F("","", eta_region.size() - 1, &eta_region[0], z_region.size() - 1, &z_region[0]) );
0390         coarse_eta_z_2D_MC[i] -> GetXaxis() -> SetTitle("#eta");
0391         coarse_eta_z_2D_MC[i] -> GetYaxis() -> SetTitle("Z vertex [mm]");
0392         coarse_eta_z_2D_MC[i] -> GetXaxis() -> SetNdivisions(505);
0393
0394         coarse_eta_z_2D_fulleta_MC.push_back( new TH2F("","", eta_region.size() - 1, &eta_region[0], z_region.size() - 1, &z_region[0]) );
0395         coarse_eta_z_2D_fulleta_MC[i] -> GetXaxis() -> SetTitle("#eta");
0396         coarse_eta_z_2D_fulleta_MC[i] -> GetYaxis() -> SetTitle("Z vertex [mm]");
0397         coarse_eta_z_2D_fulleta_MC[i] -> GetXaxis() -> SetNdivisions(505);
0398
0399         coarse_Reco_SignalNTracklet_Single_eta_z_2D.push_back( new TH2F("","", eta_region.size() - 1, &eta_region[0], z_region.size() - 1, &z_region[0]) );
0400         coarse_Reco_SignalNTracklet_Single_eta_z_2D[i] -> GetXaxis() -> SetTitle("#eta");
0401         coarse_Reco_SignalNTracklet_Single_eta_z_2D[i] -> GetYaxis() -> SetTitle("Z vertex [mm]");
0402         coarse_Reco_SignalNTracklet_Single_eta_z_2D[i] -> GetXaxis() -> SetNdivisions(505);
0403
0404         coarse_Reco_SignalNTracklet_Multi_eta_z_2D.push_back( new TH2F("","", eta_region.size() - 1, &eta_region[0], z_region.size() - 1, &z_region[0]) );
0405         coarse_Reco_SignalNTracklet_Multi_eta_z_2D[i] -> GetXaxis() -> SetTitle("#eta");
0406         coarse_Reco_SignalNTracklet_Multi_eta_z_2D[i] -> GetYaxis() -> SetTitle("Z vertex [mm]");
0407         coarse_Reco_SignalNTracklet_Multi_eta_z_2D[i] -> GetXaxis() -> SetNdivisions(505);
0408
0409         track_cluster_ratio_1D.push_back(new TH1F("","track_cluster_ratio_1D",200,0,15));
0410         track_cluster_ratio_1D[i] -> GetXaxis() -> SetTitle("NClus / NTracks");
0411         track_cluster_ratio_1D[i] -> GetYaxis() -> SetTitle("Entry");
0412         track_cluster_ratio_1D[i] -> GetXaxis() -> SetNdivisions(505);
0413
0414         track_cluster_ratio_1D_MC.push_back(new TH1F("","track_cluster_ratio_1D_MC",200,0,15));
0415         track_cluster_ratio_1D_MC[i] -> GetXaxis() -> SetTitle("NClus / NTrackMC");
0416         track_cluster_ratio_1D_MC[i] -> GetYaxis() -> SetTitle("Entry");
0417         track_cluster_ratio_1D_MC[i] -> GetXaxis() -> SetNdivisions(505);
0418
0419         dNdeta_1D.push_back(new TH1F("","",Eta_NBin, Eta_Min, Eta_Max));
0420         dNdeta_1D[i] -> SetMarkerStyle(20);
0421         dNdeta_1D[i] -> SetMarkerSize(0.8);
0422         dNdeta_1D[i] -> SetMarkerColor(TColor::GetColor("#1A3947"));
0423         dNdeta_1D[i] -> SetLineColor(TColor::GetColor("#1A3947"));
0424         dNdeta_1D[i] -> SetLineWidth(2);
0425         dNdeta_1D[i] -> GetYaxis() -> SetTitle("Entry");
0426         dNdeta_1D[i] -> GetXaxis() -> SetTitle("N tracklet #eta");
0427         // dNdeta_1D[i] -> GetYaxis() -> SetRangeUser(0,50);
0428         dNdeta_1D[i] -> SetTitleSize(0.06, "X");
0429         dNdeta_1D[i] -> SetTitleSize(0.06, "Y");
0430         dNdeta_1D[i] -> GetXaxis() -> SetTitleOffset (0.71);
0431         dNdeta_1D[i] -> GetYaxis() -> SetTitleOffset (1.1);
0432         dNdeta_1D[i] -> GetXaxis() -> CenterTitle(true);
0433         dNdeta_1D[i] -> GetYaxis() -> CenterTitle(true);
0434
0435         dNdeta_1D_MC.push_back(new TH1F("","",Eta_NBin, Eta_Min, Eta_Max));
0436         dNdeta_1D_MC[i] -> SetMarkerStyle(20);
0437         dNdeta_1D_MC[i] -> SetMarkerSize(0.8);
0438         dNdeta_1D_MC[i] -> SetMarkerColor(TColor::GetColor("#1A3947"));
0439         dNdeta_1D_MC[i] -> SetLineColor(TColor::GetColor("#1A3947"));
0440         dNdeta_1D_MC[i] -> SetLineWidth(2);
0441         dNdeta_1D_MC[i] -> GetYaxis() -> SetTitle("Entry");
0442         dNdeta_1D_MC[i] -> GetXaxis() -> SetTitle("N track #eta");
0443         // dNdeta_1D_MC[i] -> GetYaxis() -> SetRangeUser(0,50);
0444         dNdeta_1D_MC[i] -> SetTitleSize(0.06, "X");
0445         dNdeta_1D_MC[i] -> SetTitleSize(0.06, "Y");
0446         dNdeta_1D_MC[i] -> GetXaxis() -> SetTitleOffset (0.71);
0447         dNdeta_1D_MC[i] -> GetYaxis() -> SetTitleOffset (1.1);
0448         dNdeta_1D_MC[i] -> GetXaxis() -> CenterTitle(true);
0449         dNdeta_1D_MC[i] -> GetYaxis() -> CenterTitle(true);
0450
0451         dNdeta_1D_MC_edge_eta_cut.push_back(new TH1F("","",Eta_NBin, Eta_Min, Eta_Max));
0452         dNdeta_1D_MC_edge_eta_cut[i] -> SetMarkerStyle(20);
0453         dNdeta_1D_MC_edge_eta_cut[i] -> SetMarkerSize(0.8);
0454         dNdeta_1D_MC_edge_eta_cut[i] -> SetMarkerColor(TColor::GetColor("#1A3947"));
0455         dNdeta_1D_MC_edge_eta_cut[i] -> SetLineColor(TColor::GetColor("#1A3947"));
0456         dNdeta_1D_MC_edge_eta_cut[i] -> SetLineWidth(2);
0457         dNdeta_1D_MC_edge_eta_cut[i] -> GetYaxis() -> SetTitle("Entry");
0458         dNdeta_1D_MC_edge_eta_cut[i] -> GetXaxis() -> SetTitle("N track #eta");
0459         // dNdeta_1D_MC_edge_eta_cut[i] -> GetYaxis() -> SetRangeUser(0,50);
0460         dNdeta_1D_MC_edge_eta_cut[i] -> SetTitleSize(0.06, "X");
0461         dNdeta_1D_MC_edge_eta_cut[i] -> SetTitleSize(0.06, "Y");
0462         dNdeta_1D_MC_edge_eta_cut[i] -> GetXaxis() -> SetTitleOffset (0.71);
0463         dNdeta_1D_MC_edge_eta_cut[i] -> GetYaxis() -> SetTitleOffset (1.1);
0464         dNdeta_1D_MC_edge_eta_cut[i] -> GetXaxis() -> CenterTitle(true);
0465         dNdeta_1D_MC_edge_eta_cut[i] -> GetYaxis() -> CenterTitle(true);
0466
0467         final_dNdeta_1D.push_back(new TH1F("","",eta_region.size() - 1, &eta_region[0]));
0468         final_dNdeta_1D[i] -> SetMarkerStyle(20);
0469         final_dNdeta_1D[i] -> SetMarkerSize(0.8);
0470         final_dNdeta_1D[i] -> SetMarkerColor(TColor::GetColor("#1A3947"));
0471         final_dNdeta_1D[i] -> SetLineColor(TColor::GetColor("#1A3947"));
0472         final_dNdeta_1D[i] -> SetLineWidth(2);
0473         final_dNdeta_1D[i] -> GetYaxis() -> SetTitle("dN_{ch}/d#eta");
0474         final_dNdeta_1D[i] -> GetXaxis() -> SetTitle("#eta");
0475         // final_dNdeta_1D[i] -> GetYaxis() -> SetRangeUser(0,50);
0476         final_dNdeta_1D[i] -> SetTitleSize(0.06, "X");
0477         final_dNdeta_1D[i] -> SetTitleSize(0.06, "Y");
0478         final_dNdeta_1D[i] -> GetXaxis() -> SetTitleOffset (0.71);
0479         final_dNdeta_1D[i] -> GetYaxis() -> SetTitleOffset (1.1);
0480         final_dNdeta_1D[i] -> GetXaxis() -> CenterTitle(true);
0481         final_dNdeta_1D[i] -> GetYaxis() -> CenterTitle(true);
0482
0483         final_dNdeta_multi_1D.push_back(new TH1F("","",eta_region.size() - 1, &eta_region[0]));
0484         final_dNdeta_multi_1D[i] -> SetMarkerStyle(20);
0485         final_dNdeta_multi_1D[i] -> SetMarkerSize(0.8);
0486         final_dNdeta_multi_1D[i] -> SetMarkerColor(TColor::GetColor("#c48045"));
0487         final_dNdeta_multi_1D[i] -> SetLineColor(TColor::GetColor("#c48045"));
0488         final_dNdeta_multi_1D[i] -> SetLineWidth(2);
0489         final_dNdeta_multi_1D[i] -> GetYaxis() -> SetTitle("dN_{ch}/d#eta");
0490         final_dNdeta_multi_1D[i] -> GetXaxis() -> SetTitle("#eta");
0491         // final_dNdeta_multi_1D[i] -> GetYaxis() -> SetRangeUser(0,50);
0492         final_dNdeta_multi_1D[i] -> SetTitleSize(0.06, "X");
0493         final_dNdeta_multi_1D[i] -> SetTitleSize(0.06, "Y");
0494         final_dNdeta_multi_1D[i] -> GetXaxis() -> SetTitleOffset (0.71);
0495         final_dNdeta_multi_1D[i] -> GetYaxis() -> SetTitleOffset (1.1);
0496         final_dNdeta_multi_1D[i] -> GetXaxis() -> CenterTitle(true);
0497         final_dNdeta_multi_1D[i] -> GetYaxis() -> CenterTitle(true);
0498
0499         final_dNdeta_1D_MC.push_back(new TH1F("","",eta_region.size() - 1, &eta_region[0]));
0500         final_dNdeta_1D_MC[i] -> SetMarkerStyle(20);
0501         final_dNdeta_1D_MC[i] -> SetMarkerSize(0.8);
0502         final_dNdeta_1D_MC[i] -> SetMarkerColor(TColor::GetColor("#1A3947"));
0503         final_dNdeta_1D_MC[i] -> SetLineColor(TColor::GetColor("#1A3947"));
0504         final_dNdeta_1D_MC[i] -> SetLineWidth(2);
0505         final_dNdeta_1D_MC[i] -> GetYaxis() -> SetTitle("dN_{ch}/d#eta");
0506         final_dNdeta_1D_MC[i] -> GetXaxis() -> SetTitle("#eta");
0507         // final_dNdeta_1D_MC[i] -> GetYaxis() -> SetRangeUser(0,50);
0508         final_dNdeta_1D_MC[i] -> SetTitleSize(0.06, "X");
0509         final_dNdeta_1D_MC[i] -> SetTitleSize(0.06, "Y");
0510         final_dNdeta_1D_MC[i] -> GetXaxis() -> SetTitleOffset (0.71);
0511         final_dNdeta_1D_MC[i] -> GetYaxis() -> SetTitleOffset (1.1);
0512         final_dNdeta_1D_MC[i] -> GetXaxis() -> CenterTitle(true);
0513         final_dNdeta_1D_MC[i] -> GetYaxis() -> CenterTitle(true);
0514
0515         track_eta_phi_2D.push_back(new TH2F("","track_eta_phi_2D", 200, -40, 400, Eta_NBin, Eta_Min, Eta_Max));
0516         track_eta_phi_2D[i] -> GetXaxis() -> SetTitle("tracklet #phi");
0517         track_eta_phi_2D[i] -> GetYaxis() -> SetTitle("tracklet #eta");
0518         track_eta_phi_2D[i] -> GetXaxis() -> SetNdivisions(505);
0519
0520         track_eta_z_2D.push_back(new TH2F("","track_eta_z_2D", Eta_NBin, Eta_Min, Eta_Max, 300, -420, 20));
0521         track_eta_z_2D[i] -> GetXaxis() -> SetTitle("tracklet #eta");
0522         track_eta_z_2D[i] -> GetYaxis() -> SetTitle("z vertex [mm]");
0523         track_eta_z_2D[i] -> GetXaxis() -> SetNdivisions(505);
0524
0525         track_eta_z_2D_MC.push_back(new TH2F("","track_eta_z_2D_MC", Eta_NBin, Eta_Min, Eta_Max, 300, -420, 20));
0526         track_eta_z_2D_MC[i] -> GetXaxis() -> SetTitle("track #eta");
0527         track_eta_z_2D_MC[i] -> GetYaxis() -> SetTitle("z vertex [mm]");
0528         track_eta_z_2D_MC[i] -> GetXaxis() -> SetNdivisions(505);
0529
0530         track_delta_eta_1D.push_back(new TH1F("","track_delta_eta_1D", 200, -1.5, 1.5));
0531         track_delta_eta_1D[i] -> GetXaxis() -> SetTitle("Clus #Delta#eta");
0532         track_delta_eta_1D[i] -> GetYaxis() -> SetTitle("Entry");
0533         track_delta_eta_1D[i] -> GetXaxis() -> SetNdivisions(505);
0534
0535         track_delta_eta_1D_post.push_back(new TH1F("","track_delta_eta_1D_ potrack_delta_eta_1D_post", 200, -1.5, 1.5));
0536         track_delta_eta_1D_post[i] -> GetXaxis() -> SetTitle("Clus #Delta#eta");
0537         track_delta_eta_1D_post[i] -> GetYaxis() -> SetTitle("Entry");
0538         track_delta_eta_1D_post[i] -> GetXaxis() -> SetNdivisions(505);
0539
0540         track_delta_phi_1D.push_back(new TH1F("","track_delta_phi_1D", 200, -1.5, 1.5));
0541         track_delta_phi_1D[i] -> GetXaxis() -> SetTitle("Clus #Delta#phi [degree]");
0542         track_delta_phi_1D[i] -> GetYaxis() -> SetTitle("Entry");
0543         track_delta_phi_1D[i] -> GetXaxis() -> SetNdivisions(505);
0544
0545         track_DCA_distance.push_back(new TH1F("","track_DCA_distance", 200, -5, 5));
0546         track_DCA_distance[i] -> GetXaxis() -> SetTitle("DCA distance [mm]");
0547         track_DCA_distance[i] -> GetYaxis() -> SetTitle("Entry");
0548         track_DCA_distance[i] -> GetXaxis() -> SetNdivisions(505);
0549
0550         track_phi_DCA_2D.push_back(new TH2F("","track_phi_DCA_2D", 200, -1.5, 1.5, 200, -3, 3));
0551         track_phi_DCA_2D[i] -> GetXaxis() -> SetTitle("Clus #Delta#phi [degree]");
0552         track_phi_DCA_2D[i] -> GetYaxis() -> SetTitle("DCA distance [mm]");
0553         track_phi_DCA_2D[i] -> GetXaxis() -> SetNdivisions(505);
0554
0555         track_ratio_1D.push_back(new TH1F("","track_ratio_1D", 200, 0, 1.5));
0556         track_ratio_1D[i] -> GetXaxis() -> SetTitle("N reco track / N true track");
0557         track_ratio_1D[i] -> GetYaxis() -> SetTitle("Entry");
0558         track_ratio_1D[i] -> GetXaxis() -> SetNdivisions(505);
0559
0560         track_DeltaPhi_eta_2D.push_back(new TH2F("","track_DeltaPhi_eta_2D", 200, -1.5, 1.5, 200, -4, 4));
0561         track_DeltaPhi_eta_2D[i] -> GetXaxis() -> SetTitle("Clus #Delta#phi [degree]");
0562         track_DeltaPhi_eta_2D[i] -> GetYaxis() -> SetTitle("Tracklet #eta");
0563         track_DeltaPhi_eta_2D[i] -> GetXaxis() -> SetNdivisions(505);
0564
0565         track_DeltaPhi_DeltaEta_2D.push_back(new TH2F("","track_DeltaPhi_DeltaEta_2D", 200, -1.5, 1.5, 200, -1, 1));
0566         track_DeltaPhi_DeltaEta_2D[i] -> GetXaxis() -> SetTitle("Clus #Delta#phi [degree]");
0567         track_DeltaPhi_DeltaEta_2D[i] -> GetYaxis() -> SetTitle("#Delta#eta");
0568         track_DeltaPhi_DeltaEta_2D[i] -> GetXaxis() -> SetNdivisions(505);
0569
0570         final_track_delta_phi_1D[i].clear();
0571         final_track_multi_delta_phi_1D[i].clear();
0572         // for (int eta_i = 0; eta_i < eta_region.size() - 1; eta_i++)
0573         // note : now we try to do something inclusive, which means we open the z selection, try to have everything, and then we decide afterwards
0574         for (int eta_i = 0; eta_i < (coarse_eta_z_map -> GetNbinsX() * coarse_eta_z_map -> GetNbinsY()); eta_i++)
0575         {
0576             final_track_delta_phi_1D[i].push_back(new TH1F("","final_track_delta_phi_1D", 100, -3, 3));
0577             final_track_delta_phi_1D[i][eta_i] -> GetXaxis() -> SetTitle("Clus #Delta#phi [degree]");
0578             final_track_delta_phi_1D[i][eta_i] -> GetYaxis() -> SetTitle("Entry");
0579             final_track_delta_phi_1D[i][eta_i] -> GetXaxis() -> SetNdivisions(505);
0580
0581             final_track_multi_delta_phi_1D[i].push_back(new TH1F("","final_track_multi_delta_phi_1D", 100, -3, 3));
0582             final_track_multi_delta_phi_1D[i][eta_i] -> GetXaxis() -> SetTitle("Clus #Delta#phi [degree]");
0583             final_track_multi_delta_phi_1D[i][eta_i] -> GetYaxis() -> SetTitle("Entry");
0584             final_track_multi_delta_phi_1D[i][eta_i] -> GetXaxis() -> SetNdivisions(505);
0585         }
0586     }
0587
0588     mega_track_contamination_2D = new TH2F("","", 200,0,8500, 200, 0, 5);
0589     mega_track_contamination_2D -> GetXaxis() -> SetTitle("INTT NClus");
0590     mega_track_contamination_2D -> GetYaxis() -> SetTitle("Mega Track contamination ratio [%]");
0591     mega_track_contamination_2D -> GetXaxis() -> SetNdivisions(505);
0592
0593     mega_track_finding_ratio_2D = new TH2F("","", 200,0,8500, 200, 0, 2);
0594     mega_track_finding_ratio_2D -> GetXaxis() -> SetTitle("INTT NClus");
0595     mega_track_finding_ratio_2D -> GetYaxis() -> SetTitle("Mega Track finding ratio");
0596     mega_track_finding_ratio_2D -> GetXaxis() -> SetNdivisions(505);
0597
0598     MBin_Z_evt_hist_2D = new TH2F("","", centrality_region.size(), 0, centrality_region.size(), z_region.size() - 1, &z_region[0]);
0599     MBin_Z_evt_hist_2D -> GetXaxis() -> SetTitle("Centrality bin");
0600     MBin_Z_evt_hist_2D -> GetYaxis() -> SetTitle("Z vertex [mm]");
0601     MBin_Z_evt_hist_2D -> GetXaxis() -> SetNdivisions(505);
0602
0603     MBin_Z_evt_hist_2D_MC = new TH2F("","", centrality_region.size(), 0, centrality_region.size(), z_region.size() - 1, &z_region[0]);
0604     MBin_Z_evt_hist_2D_MC -> GetXaxis() -> SetTitle("Centrality bin");
0605     MBin_Z_evt_hist_2D_MC -> GetYaxis() -> SetTitle("Z vertex [mm]");
0606     MBin_Z_evt_hist_2D_MC -> GetXaxis() -> SetNdivisions(505);
0607
0608
0609     track_cluster_ratio_multiplicity_2D = new TH2F("","track_cluster_ratio_multiplicity_2D", 200, 0, 9000, 200, 0, 15);
0610     track_cluster_ratio_multiplicity_2D -> GetXaxis() -> SetTitle("NClus");
0611     track_cluster_ratio_multiplicity_2D -> GetYaxis() -> SetTitle("NClus / NTracks");
0612     track_cluster_ratio_multiplicity_2D -> GetXaxis() -> SetNdivisions(505);
0613
0614     track_cluster_ratio_multiplicity_2D_MC = new TH2F("","track_cluster_ratio_multiplicity_2D_MC", 200, 0, 9000, 200, 0, 15);
0615     track_cluster_ratio_multiplicity_2D_MC -> GetXaxis() -> SetTitle("NClus");
0616     track_cluster_ratio_multiplicity_2D_MC -> GetYaxis() -> SetTitle("NClus / NTracks");
0617     track_cluster_ratio_multiplicity_2D_MC -> GetXaxis() -> SetNdivisions(505);
0618
0619     track_correlation_2D = new TH2F("","track_correlation_2D", 200, 0, 3000, 200, 0, 3000);
0620     track_correlation_2D -> GetXaxis() -> SetTitle("N true track");
0621     track_correlation_2D -> GetYaxis() -> SetTitle("N reco track");
0622     track_correlation_2D -> GetXaxis() -> SetNdivisions(505);
0623
0624     track_ratio_2D = new TH2F("","track_ratio_2D", 200, 0, 3000, 200, 0, 2);
0625     track_ratio_2D -> GetXaxis() -> SetTitle("N true track");
0626     track_ratio_2D -> GetYaxis() -> SetTitle("N reco track/N true track");
0627     track_ratio_2D -> GetXaxis() -> SetNdivisions(505);
0628
0629     check_inner_layer_clu_phi_1D = new TH1F("","check_inner_layer_clu_phi_1D", 440, -40, 400);
0630     check_inner_layer_clu_phi_1D -> GetXaxis() -> SetTitle("Clus #phi [degree]");
0631     check_inner_layer_clu_phi_1D -> GetYaxis() -> SetTitle("Entry");
0632     // check_inner_layer_clu_phi_1D -> GetXaxis() -> SetNdivisions(505);
0633
0634     check_outer_layer_clu_phi_1D = new TH1F("","check_outer_layer_clu_phi_1D", 440, -40, 400);
0635     check_outer_layer_clu_phi_1D -> GetXaxis() -> SetTitle("Clus #phi [degree]");
0636     check_outer_layer_clu_phi_1D -> GetYaxis() -> SetTitle("Entry");
0637     // check_outer_layer_clu_phi_1D -> GetXaxis() -> SetNdivisions(505);
0638
0639     reco_eta_correlation_2D = new TH2F("","reco_eta_correlation_2D", 200, -3, 3, 200, -3, 3);
0640     reco_eta_correlation_2D -> GetXaxis() -> SetTitle("Reco 3P #eta");
0641     reco_eta_correlation_2D -> GetYaxis() -> SetTitle("Reco 2P #eta");
0642     reco_eta_correlation_2D -> GetXaxis() -> SetNdivisions(505);
0643
0644     reco_eta_diff_reco3P_2D = new TH2F("","reco_eta_diff_reco3P_2D", 200, -3, 3, 200, -0.35, 0.35);
0645     reco_eta_diff_reco3P_2D -> GetXaxis() -> SetTitle("Reco 3P #eta");
0646     reco_eta_diff_reco3P_2D -> GetYaxis() -> SetTitle("Reco #eta diff");
0647     reco_eta_diff_reco3P_2D -> GetXaxis() -> SetNdivisions(505);
0648
0649     reco_eta_diff_1D = new TH1F("","reco_eta_diff_1D", 200, -0.35, 0.35);
0650     reco_eta_diff_1D -> GetXaxis() -> SetTitle("Reco #eta diff");
0651     reco_eta_diff_1D -> GetYaxis() -> SetTitle("Entry");
0652     reco_eta_diff_1D -> GetXaxis() -> SetNdivisions(505);
0653
0654     clu_used_centrality_2D = new TH2F("","clu_used_centrality_2D", 200, 0, 8500, 10, 0, 10);
0655     clu_used_centrality_2D -> GetXaxis() -> SetTitle("NClus");
0656     clu_used_centrality_2D -> GetYaxis() -> SetTitle("N used each clu");
0657     clu_used_centrality_2D -> GetXaxis() -> SetNdivisions(505);
0658
0659     exclusive_NClus_inner        = new TH1F("","exclusive_NClus_inner;NClus inner barrel;Entry",50, 0, 5000);
0660     exclusive_NClus_outer        = new TH1F("","exclusive_NClus_outer;NClus outer barrel;Entry",50,0,5000);
0661     exclusive_NClus_sum          = new TH1F("","exclusive_NClus_sum;NClus total;Entry",50,0,9000);
0662
0663     exclusive_cluster_inner_eta  = new TH1F("","exclusive_cluster_inner_eta;Cluster inner #eta;Entry",50,-3.5,3.5);
0664     exclusive_cluster_inner_phi  = new TH1F("","exclusive_cluster_inner_phi;Cluster inner #phi [radian];Entry",50,-3.5,3.5);
0665     exclusive_cluster_outer_eta  = new TH1F("","exclusive_cluster_outer_eta;Cluster outer #eta;Entry",50,-3.5,3.5);
0666     exclusive_cluster_outer_phi  = new TH1F("","exclusive_cluster_outer_phi;Cluster outer #phi [radian];Entry",50,-3.5,3.5);
0667     exclusive_cluster_all_eta    = new TH1F("","exclusive_cluster_all_eta;Cluster #eta;Entry",50,-3.5,3.5);
0668     exclusive_cluster_all_phi    = new TH1F("","exclusive_cluster_all_phi;Cluster #phi [radian];Entry",50,-3.5,3.5);
0669
0670     exclusive_cluster_inner_adc  = new TH1F("","exclusive_cluster_inner_adc;Cluster inner adc;Entry",50,0,400);
0671     exclusive_cluster_outer_adc  = new TH1F("","exclusive_cluster_outer_adc;Cluster outer adc;Entry",50,0,400);
0672
0673     exclusive_cluster_inner_eta_adc_2D = new TH2F("","exclusive_cluster_inner_eta_adc_2D;Cluster inner #eta; Cluster ADC",
0674         2000, exclusive_cluster_inner_eta -> GetXaxis() -> GetXmin(), exclusive_cluster_inner_eta -> GetXaxis() -> GetXmax(),
0675         exclusive_cluster_inner_adc -> GetNbinsX(), exclusive_cluster_inner_adc -> GetXaxis() -> GetXmin(), exclusive_cluster_inner_adc -> GetXaxis() -> GetXmax()
0676     );
0677
0678     exclusive_cluster_outer_eta_adc_2D = new TH2F("","exclusive_cluster_outer_eta_adc_2D;Cluster outer #eta; Cluster ADC",
0679         2000, exclusive_cluster_outer_eta -> GetXaxis() -> GetXmin(), exclusive_cluster_outer_eta -> GetXaxis() -> GetXmax(),
0680         exclusive_cluster_outer_adc -> GetNbinsX(), exclusive_cluster_outer_adc -> GetXaxis() -> GetXmin(), exclusive_cluster_outer_adc -> GetXaxis() -> GetXmax()
0681     );
0682
0683     exclusive_cluster_inner_eta_phi_2D = new TH2F(
0684         "",
0685         "exclusive_cluster_inner_eta_phi_2D;Cluster inner #eta; Cluster inner #phi",
0686         2000, exclusive_cluster_inner_eta -> GetXaxis() -> GetXmin(), exclusive_cluster_inner_eta -> GetXaxis() -> GetXmax(),
0687         2000, exclusive_cluster_inner_phi -> GetXaxis() -> GetXmin(), exclusive_cluster_inner_phi -> GetXaxis() -> GetXmax()
0688     );
0689     exclusive_cluster_outer_eta_phi_2D = new TH2F(
0690         "",
0691         "exclusive_cluster_outer_eta_phi_2D;Cluster outer #eta; Cluster outer #phi",
0692         2000, exclusive_cluster_outer_eta -> GetXaxis() -> GetXmin(), exclusive_cluster_outer_eta -> GetXaxis() -> GetXmax(),
0693         2000, exclusive_cluster_outer_phi -> GetXaxis() -> GetXmin(), exclusive_cluster_outer_phi -> GetXaxis() -> GetXmax()
0694     );
0695
0696     exclusive_tight_tracklet_eta = new TH1F("","exclusive_tight_tracklet_eta;Tracklet #eta;Entry",50,-3.5,3.5);
0697     exclusive_tight_tracklet_phi = new TH1F("","exclusive_tight_tracklet_phi;Tracklet #phi [radian];Entry",50,-3.5,3.5);
0698     exclusive_loose_tracklet_eta = new TH1F("","exclusive_loose_tracklet_eta;Tracklet #eta;Entry",50,-3.5,3.5);
0699     exclusive_loose_tracklet_phi = new TH1F("","exclusive_loose_tracklet_phi;Tracklet #phi [radian];Entry",50,-3.5,3.5);
0700
0701     eta_region_hist = new TH1F("","", eta_region.size() - 1, &eta_region[0]);
0702
0703     cout<<" running ? in INTTEta, InitHist"<<endl;
0704
0705 }
0706
0707 void INTTEta::InitGraph()
0708 {
0709     track_gr = new TGraphErrors();
0710
0711     evt_reco_track_gr_All = new TGraph();
0712     evt_reco_track_gr_All -> SetMarkerStyle(5);
0713     evt_reco_track_gr_All -> SetMarkerSize(1);
0714     evt_reco_track_gr_All -> SetMarkerColor(4);
0715
0716     evt_reco_track_gr_PhiLoose = new TGraph();
0717     evt_reco_track_gr_PhiLoose -> SetMarkerStyle(5);
0718     evt_reco_track_gr_PhiLoose -> SetMarkerSize(1);
0719     evt_reco_track_gr_PhiLoose -> SetMarkerColor(4);
0720
0721     evt_reco_track_gr_Z = new TGraph();
0722     evt_reco_track_gr_Z -> SetMarkerStyle(5);
0723     evt_reco_track_gr_Z -> SetMarkerSize(1);
0724     evt_reco_track_gr_Z -> SetMarkerColor(4);
0725
0726     evt_reco_track_gr_ZDCA = new TGraph();
0727     evt_reco_track_gr_ZDCA -> SetMarkerStyle(5);
0728     evt_reco_track_gr_ZDCA -> SetMarkerSize(1);
0729     evt_reco_track_gr_ZDCA -> SetMarkerColor(4);
0730
0731     evt_reco_track_gr_ZDCAPhi = new TGraph();
0732     evt_reco_track_gr_ZDCAPhi -> SetMarkerStyle(5);
0733     evt_reco_track_gr_ZDCAPhi -> SetMarkerSize(1);
0734     evt_reco_track_gr_ZDCAPhi -> SetMarkerColor(4);
0735
0736     evt_true_track_gr = new TGraph();
0737     evt_true_track_gr -> SetMarkerStyle(20);
0738     evt_true_track_gr -> SetMarkerSize(1);
0739     // evt_true_track_gr -> SetMarkerColor(2);
0740     evt_true_track_gr -> SetMarkerColorAlpha(2,0.5);
0741
0742
0743     cout<<" running ? "<<endl;
0744     return;
0745 }
0746
0747 void INTTEta::InitTree()
0748 {
0749     out_file = new TFile(Form("%s/INTT_final_hist_info.root",out_folder_directory.c_str()),"RECREATE");
0750     tree_out =  new TTree ("tree_eta", "Tracklet eta info.");
0751
0752     tree_out -> Branch("eID",&out_eID);
0753     tree_out -> Branch("Evt_centrality_bin",&out_evt_centrality_bin);
0754     tree_out -> Branch("Evt_zvtx", &out_evt_zvtx);
0755     tree_out -> Branch("True_zvtx", &out_true_zvtx);
0756     tree_out -> Branch("NTrueTrack", &out_NTrueTrack);
0757     tree_out -> Branch("NClus", &out_total_NClus);
0758     tree_out -> Branch("RecoTrack_eta_d", &out_recotrack_eta_d);
0759     tree_out -> Branch("RecoTrack_phi_d", &out_recotrack_phi_d);
0760     tree_out -> Branch("TrueTrack_eta_d", &out_truetrack_eta_d);
0761     tree_out -> Branch("TrueTrack_phi_d", &out_truetrack_phi_d);
0762     // tree_out -> Branch("RecoTrack_eta_i", &out_track_eta_i);
0763     // tree_out -> Branch("Track_delta_phi_d", &out_track_delta_phi_d);
0764     tree_out -> Branch("N2Clu_track", &out_N2Clu_track);
0765     tree_out -> Branch("N3Clu_track", &out_N3Clu_track);
0766     tree_out -> Branch("N4Clu_track", &out_N4Clu_track);
0767
0768     return;
0769 }
0770
0771 void INTTEta::print_evt_plot(int event_i, int NTrueTrack, int innerNClu, int outerNClu)
0772 {
0773     c1 -> Clear();
0774     c1 -> Print( Form("%s/evt_track_eID_%i.pdf(", out_folder_directory_evt.c_str(), event_i) );
0775     c1 -> cd();
0776
0777     evt_reco_track_gr_All -> GetXaxis() -> SetTitle("#eta");
0778     evt_reco_track_gr_All -> GetYaxis() -> SetTitle("#phi");
0779     evt_reco_track_gr_All -> GetXaxis() -> SetLimits(-3.5, 3.5);
0780     evt_reco_track_gr_All -> GetYaxis() -> SetRangeUser(-10, 420);
0781     evt_reco_track_gr_All -> GetXaxis() -> SetNdivisions(505);
0782     evt_reco_track_gr_All -> Draw("ap");
0783     evt_true_track_gr -> Draw("p same");
0784     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0785     draw_text -> DrawLatex(0.21, 0.90, Form("NTrueTrack: %i, innerNClu: %i, outerNClu: %i", NTrueTrack, innerNClu, outerNClu));
0786     draw_text -> DrawLatex(0.21, 0.85, Form("NReco_tracklet: %i", evt_reco_track_gr_All->GetN()));
0787     coord_line -> DrawLine(-3.5, 0, 3.5, 0);
0788     coord_line -> DrawLine(-3.5, 360, 3.5, 360);
0789     c1 -> Print( Form("%s/evt_track_eID_%i.pdf", out_folder_directory_evt.c_str(), event_i) );
0790     c1 -> Clear();
0791
0792     evt_reco_track_gr_PhiLoose -> GetXaxis() -> SetTitle("#eta");
0793     evt_reco_track_gr_PhiLoose -> GetYaxis() -> SetTitle("#phi");
0794     evt_reco_track_gr_PhiLoose -> GetXaxis() -> SetLimits(-3.5, 3.5);
0795     evt_reco_track_gr_PhiLoose -> GetYaxis() -> SetRangeUser(-10, 420);
0796     evt_reco_track_gr_PhiLoose -> GetXaxis() -> SetNdivisions(505);
0797     evt_reco_track_gr_PhiLoose -> Draw("ap");
0798     evt_true_track_gr -> Draw("p same");
0799     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0800     draw_text -> DrawLatex(0.21, 0.90, Form("NTrueTrack: %i, innerNClu: %i, outerNClu: %i", NTrueTrack, innerNClu, outerNClu));
0801     draw_text -> DrawLatex(0.21, 0.85, Form("NReco_tracklet: %i", evt_reco_track_gr_PhiLoose->GetN()));
0802     coord_line -> DrawLine(-3.5, 0, 3.5, 0);
0803     coord_line -> DrawLine(-3.5, 360, 3.5, 360);
0804     c1 -> Print( Form("%s/evt_track_eID_%i.pdf", out_folder_directory_evt.c_str(), event_i) );
0805     c1 -> Clear();
0806
0807     evt_reco_track_gr_Z -> GetXaxis() -> SetTitle("#eta");
0808     evt_reco_track_gr_Z -> GetYaxis() -> SetTitle("#phi");
0809     evt_reco_track_gr_Z -> GetXaxis() -> SetLimits(-3.5, 3.5);
0810     evt_reco_track_gr_Z -> GetYaxis() -> SetRangeUser(-10, 420);
0811     evt_reco_track_gr_Z -> GetXaxis() -> SetNdivisions(505);
0812     evt_reco_track_gr_Z -> Draw("ap");
0813     evt_true_track_gr -> Draw("p same");
0814     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0815     draw_text -> DrawLatex(0.21, 0.90, Form("NTrueTrack: %i, innerNClu: %i, outerNClu: %i", NTrueTrack, innerNClu, outerNClu));
0816     draw_text -> DrawLatex(0.21, 0.85, Form("NReco_tracklet: %i", evt_reco_track_gr_Z->GetN()));
0817     coord_line -> DrawLine(-3.5, 0, 3.5, 0);
0818     coord_line -> DrawLine(-3.5, 360, 3.5, 360);
0819     c1 -> Print( Form("%s/evt_track_eID_%i.pdf", out_folder_directory_evt.c_str(), event_i) );
0820     c1 -> Clear();
0821
0822     evt_reco_track_gr_ZDCA -> GetXaxis() -> SetTitle("#eta");
0823     evt_reco_track_gr_ZDCA -> GetYaxis() -> SetTitle("#phi");
0824     evt_reco_track_gr_ZDCA -> GetXaxis() -> SetLimits(-3.5, 3.5);
0825     evt_reco_track_gr_ZDCA -> GetYaxis() -> SetRangeUser(-10, 420);
0826     evt_reco_track_gr_ZDCA -> GetXaxis() -> SetNdivisions(505);
0827     evt_reco_track_gr_ZDCA -> Draw("ap");
0828     evt_true_track_gr -> Draw("p same");
0829     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0830     draw_text -> DrawLatex(0.21, 0.90, Form("NTrueTrack: %i, innerNClu: %i, outerNClu: %i", NTrueTrack, innerNClu, outerNClu));
0831     draw_text -> DrawLatex(0.21, 0.85, Form("NReco_tracklet: %i", evt_reco_track_gr_ZDCA->GetN()));
0832     coord_line -> DrawLine(-3.5, 0, 3.5, 0);
0833     coord_line -> DrawLine(-3.5, 360, 3.5, 360);
0834     c1 -> Print( Form("%s/evt_track_eID_%i.pdf", out_folder_directory_evt.c_str(), event_i) );
0835     c1 -> Clear();
0836
0837     evt_reco_track_gr_ZDCAPhi -> GetXaxis() -> SetTitle("#eta");
0838     evt_reco_track_gr_ZDCAPhi -> GetYaxis() -> SetTitle("#phi");
0839     evt_reco_track_gr_ZDCAPhi -> GetXaxis() -> SetLimits(-3.5, 3.5);
0840     evt_reco_track_gr_ZDCAPhi -> GetYaxis() -> SetRangeUser(-10, 420);
0841     evt_reco_track_gr_ZDCAPhi -> GetXaxis() -> SetNdivisions(505);
0842     evt_reco_track_gr_ZDCAPhi -> Draw("ap");
0843     evt_true_track_gr -> Draw("p same");
0844     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
0845     draw_text -> DrawLatex(0.21, 0.90, Form("NTrueTrack: %i, innerNClu: %i, outerNClu: %i", NTrueTrack, innerNClu, outerNClu));
0846     draw_text -> DrawLatex(0.21, 0.85, Form("NReco_tracklet: %i", evt_reco_track_gr_ZDCAPhi->GetN()));
0847     coord_line -> DrawLine(-3.5, 0, 3.5, 0);
0848     coord_line -> DrawLine(-3.5, 360, 3.5, 360);
0849     c1 -> Print( Form("%s/evt_track_eID_%i.pdf)", out_folder_directory_evt.c_str(), event_i) );
0850     c1 -> Clear();
0851 }
0852
0853 // // note : this function is for the event by event vertex calculation
0854 // // note : this function is the old method, which means there are two for loops that checks all the combinations, it takes a lot of time
0855 // void INTTEta::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, int centrality_bin, vector<vector<float>> true_track_info) // note : evt_z : {z, width}
0856 // {
0857 //     return_tag = 0;
0858
0859 //     if (event_i%1 == 0) {cout<<"In INTTEta class, running event : "<<event_i<<endl;}
0860
0861 //     inner_NClu = temp_sPH_inner_nocolumn_vec.size();
0862 //     outer_NClu = temp_sPH_outer_nocolumn_vec.size();
0863 //     total_NClus = inner_NClu + outer_NClu;
0864
0865 //     // cout<<"test_0"<<endl;
0866 //     if (total_NClus < zvtx_cal_require) {return; cout<<"return confirmation"<<endl;}
0867
0868 //     if (run_type == "MC" && NvtxMC != 1) { return; cout<<"In INTTEta class, event : "<<event_i<<" Nvtx : "<<NvtxMC<<" Nvtx more than one "<<endl;}
0869 //     if (PhiCheckTag == false)            { return; cout<<"In INTTEta class, event : "<<event_i<<" Nvtx : "<<NvtxMC<<" Not full phi has hits "<<endl;}
0870
0871 //     if (inner_NClu < 10 || outer_NClu < 10 || total_NClus > N_clu_cut || total_NClus < N_clu_cutl)
0872 //     {
0873 //         return;
0874 //         printf("In INTTEta class, event : %i, low clu continue, NClus : %lu \n", event_i, total_NClus);
0875 //     }
0876
0877 //     // todo : the z vertex range is here
0878 //     if (-220 > evt_z.first + evt_z.second || -180 < evt_z.first - evt_z.second) {return;}
0879
0880
0881 //     N_GoodEvent += 1;
0882 //     N_GoodEvent_vec[centrality_map[centrality_bin]] += 1;
0883
0884 //     // cout<<"N inner cluster : "<<inner_NClu<<" N outer cluster : "<<outer_NClu<<endl;
0885
0886 //     if (run_type == "MC")
0887 //     {
0888 //         // note : for the true track case ----------------------------------------------------------------------------------------------------------------------------------------------------------------
0889 //         double INTT_eta_acceptance_l = -0.5 * TMath::Log((sqrt(pow(-230.-TrigvtxMC[2]*10.,2)+pow(INTT_layer_R[3],2))-(-230.-TrigvtxMC[2]*10.)) / (sqrt(pow(-230.-TrigvtxMC[2]*10.,2)+pow(INTT_layer_R[3],2))+(-230.-TrigvtxMC[2]*10.))); // note : left
0890 //         double INTT_eta_acceptance_r =  -0.5 * TMath::Log((sqrt(pow(230.-TrigvtxMC[2]*10.,2)+pow(INTT_layer_R[3],2))-(230.-TrigvtxMC[2]*10.)) / (sqrt(pow(230.-TrigvtxMC[2]*10.,2)+pow(INTT_layer_R[3],2))+(230.-TrigvtxMC[2]*10.))); // note : right
0891 //         // if (event_i % 100 == 0){cout<<"z : "<<TrigvtxMC[2]*10.<<" eta : "<<INTT_eta_acceptance_l<<" "<<INTT_eta_acceptance_r<<endl;}
0892
0893 //         // cout<<"true_track_info : "<<true_track_info.size()<<endl;
0894 //         for (int track_i = 0; track_i < true_track_info.size(); track_i++)
0895 //         {
0896 //             if (true_track_info[track_i][2] == 111 || true_track_info[track_i][2] == 22 || abs(true_track_info[track_i][2]) == 2112){continue;}
0897
0898 //             if (true_track_info[track_i][0] > INTT_eta_acceptance_l && true_track_info[track_i][0] < INTT_eta_acceptance_r)
0899 //             {
0900 //                 dNdeta_1D_MC[centrality_map[centrality_bin]] -> Fill(true_track_info[track_i][0]);
0901 //                 dNdeta_1D_MC[dNdeta_1D_MC.size() - 1]        -> Fill(true_track_info[track_i][0]);
0902 //                 final_dNdeta_1D_MC[centrality_map[centrality_bin]] -> Fill(true_track_info[track_i][0]);
0903 //                 final_dNdeta_1D_MC[dNdeta_1D_MC.size() - 1]        -> Fill(true_track_info[track_i][0]);
0904 //                 track_eta_z_2D_MC[centrality_map[centrality_bin]] -> Fill(true_track_info[track_i][0], TrigvtxMC[2]*10.);
0905 //                 track_eta_z_2D_MC[track_eta_z_2D_MC.size() - 1]   -> Fill(true_track_info[track_i][0], TrigvtxMC[2]*10.);
0906
0907 //                 // cout<<"true track eta : "<<true_track_info[track_i][0]<<" phi : "<<convertTo360(true_track_info[track_i][1])<<endl;
0908 //                 // cout<<"("<<true_track_info[track_i][0]<<", "<<convertTo360(true_track_info[track_i][1])<<"), ";
0909
0910 //                 evt_true_track_gr -> SetPoint(evt_true_track_gr->GetN(),true_track_info[track_i][0], convertTo360(true_track_info[track_i][1]));
0911
0912 //                 evt_NTrack_MC += 1;
0913 //             }
0914 //         }
0915 //         // if (evt_NTrack_MC < 10) {cout<<"evt : "<<event_i<<" ---- N reco track : "<<evt_NTrack<<" N true track : "<<evt_NTrack_MC<<" ratio : "<<double(evt_NTrack) / double(evt_NTrack_MC)<<endl;}
0916 //     }
0917
0918
0919 //     // cout<<"test_1"<<endl;
0920 //     for ( int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++ )
0921 //     {
0922 //         loop_NGoodPair = 0;
0923
0924 //         for ( int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++ )
0925 //         {
0926 //             // cout<<event_i<<" test_5 "<<inner_i<<" "<<outer_i<<endl;
0927 //             // note : calculate the cluster phi after the vertex correction which can enhence the purity of the tracklet selection
0928 //             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());
0929 //             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());
0930 //             double delta_phi = Clus_InnerPhi_Offset - Clus_OuterPhi_Offset;
0931 //             double track_phi = (Clus_InnerPhi_Offset + Clus_OuterPhi_Offset)/2.;
0932 //             double inner_clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y, temp_sPH_inner_nocolumn_vec[inner_i].z});
0933 //             double outer_clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y, temp_sPH_outer_nocolumn_vec[outer_i].z});
0934 //             double delta_eta = inner_clu_eta - outer_clu_eta;
0935
0936 //             // cout<<"test_2"<<endl;
0937 //             // note : before calculating the possible z vertex range of one tracklet, the vertex correction was applied
0938 //             pair<double,double> z_range_info = Get_possible_zvtx(
0939 //                 0., // get_radius(beam_origin.first,beam_origin.second),
0940 //                 {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
0941 //                 {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
0942 //             );
0943
0944 //             if (run_type == "MC" && evt_NTrack_MC < draw_evt_cut && event_i % print_plot_ratio == 0)
0945 //             {
0946 //                 Get_eta_pair = Get_eta(
0947 //                     {get_radius(beam_origin.first,beam_origin.second), evt_z.first,evt_z.second},
0948 //                     {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},
0949 //                     {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}
0950 //                 );
0951
0952 //                 evt_reco_track_gr_All -> SetPoint(evt_reco_track_gr_All->GetN(),Get_eta_pair.second, track_phi);
0953 //             }
0954
0955
0956 //             // if ( Get_eta_pair.second > 1.4 && Get_eta_pair.second < 1.7 && track_phi > 50 && track_phi < 60 ) {
0957 //             //     cout<<" "<<endl;
0958 //             //     cout<<"reco eta : "<<Get_eta_pair.second<<" reco phi : "<<track_phi<<" delta eta : "<<delta_eta<<endl;
0959 //             //     cout<<"all pair info, inner clu pos "<<temp_sPH_inner_nocolumn_vec[inner_i].x<<" "<<temp_sPH_inner_nocolumn_vec[inner_i].y<<" "<<temp_sPH_inner_nocolumn_vec[inner_i].z<<" phi : "<<Clus_InnerPhi_Offset<<endl;
0960 //             //     cout<<"all pair info, outer clu pos "<<temp_sPH_outer_nocolumn_vec[outer_i].x<<" "<<temp_sPH_outer_nocolumn_vec[outer_i].y<<" "<<temp_sPH_outer_nocolumn_vec[outer_i].z<<" phi : "<<Clus_OuterPhi_Offset<<endl;
0961 //             //     cout<<"z vertex range : "<<z_range_info.first - z_range_info.second<<" "<<z_range_info.first + z_range_info.second<<endl;
0962 //             //     cout<<"accept z range : "<<evt_z.first - evt_z.second<<" "<<evt_z.first + evt_z.second<<endl;
0963
0964 //             // }
0965
0966 //             if (fabs(delta_phi) > 5.72) {continue;}
0967 //             if (run_type == "MC" && evt_NTrack_MC < draw_evt_cut && event_i % print_plot_ratio == 0) {evt_reco_track_gr_PhiLoose -> SetPoint(evt_reco_track_gr_PhiLoose->GetN(),Get_eta_pair.second, track_phi);}
0968
0969 //             track_delta_eta_1D[centrality_map[centrality_bin]] -> Fill( delta_eta );
0970 //             track_delta_eta_1D[track_delta_eta_1D.size() - 1]  -> Fill( delta_eta );
0971
0972 //             track_DeltaPhi_DeltaEta_2D[centrality_map[centrality_bin]]        -> Fill(delta_phi, delta_eta);
0973 //             track_DeltaPhi_DeltaEta_2D[track_DeltaPhi_DeltaEta_2D.size() - 1] -> Fill(delta_phi, delta_eta);
0974
0975 //             // cout<<event_i<<" test_6 "<<inner_i<<" "<<outer_i<<endl;
0976 //             // 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
0977 //             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;}
0978 //             // 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;
0979 //             // cout<<"test_3"<<endl;
0980
0981 //             if (run_type == "MC" && evt_NTrack_MC < draw_evt_cut && event_i % print_plot_ratio == 0) {evt_reco_track_gr_Z -> SetPoint(evt_reco_track_gr_Z->GetN(),Get_eta_pair.second, track_phi);}
0982
0983 //             // note : the reason to use the DCA calculation with sign is because that the distribution of 1D DCA will be single peak only
0984 //             double DCA_sign = calculateAngleBetweenVectors(
0985 //                 temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y,
0986 //                 temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y,
0987 //                 beam_origin.first, beam_origin.second
0988 //             );
0989
0990 //             // vector<double> DCA_info_vec = calculateDistanceAndClosestPoint(
0991 //             //     temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y,
0992 //             //     temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y,
0993 //             //     beam_origin.first, beam_origin.second
0994 //             // );
0995 //             // if (DCA_info_vec[0] != fabs(DCA_sign) && fabs( DCA_info_vec[0] - fabs(DCA_sign) ) > 0.1 ){
0996 //             //     cout<<"different DCA : "<<DCA_info_vec[0]<<" "<<DCA_sign<<" diff : "<<DCA_info_vec[0] - fabs(DCA_sign)<<endl;
0997 //             // }
0998
0999 //             track_delta_eta_1D_post[centrality_map[centrality_bin]]     -> Fill( delta_eta );
1000 //             track_delta_eta_1D_post[track_delta_eta_1D_post.size() - 1] -> Fill( delta_eta );
1001
1002 //             track_delta_phi_1D[centrality_map[centrality_bin]] -> Fill( delta_phi );
1003 //             track_delta_phi_1D[track_delta_phi_1D.size() - 1]  -> Fill( delta_phi );
1004
1005 //             track_DCA_distance[centrality_map[centrality_bin]] -> Fill( DCA_sign );
1006 //             track_DCA_distance[track_DCA_distance.size() - 1]  -> Fill( DCA_sign );
1007
1008 //             track_phi_DCA_2D[centrality_map[centrality_bin]] -> Fill( delta_phi, DCA_sign );
1009 //             track_phi_DCA_2D[track_phi_DCA_2D.size() - 1]    -> Fill( delta_phi, DCA_sign );
1010
1011 //             // cout<<" "<<endl;
1012 //             // cout<<"inner_i : "<<inner_i<<" outer_i "<<outer_i<<" true z : "<<TrigvtxMC[2]*10.<<" reco z: "<<evt_z.first<<endl;
1013 //             // cout<<"all pair info, inner clu pos "<<temp_sPH_inner_nocolumn_vec[inner_i].x<<" "<<temp_sPH_inner_nocolumn_vec[inner_i].y<<" "<<temp_sPH_inner_nocolumn_vec[inner_i].z<<" phi : "<<Clus_InnerPhi_Offset<<endl;
1014 //             // cout<<"all pair info, outer clu pos "<<temp_sPH_outer_nocolumn_vec[outer_i].x<<" "<<temp_sPH_outer_nocolumn_vec[outer_i].y<<" "<<temp_sPH_outer_nocolumn_vec[outer_i].z<<" phi : "<<Clus_OuterPhi_Offset<<endl;
1015 //             // if (fabs(delta_phi) < 4) {cout<<"("<<Get_eta_pair.second<<", "<<track_phi<<"), ";}
1016
1017 //             // if (DCA_cut.first < DCA_sign && DCA_sign < DCA_cut.second)
1018 //             if (true)
1019 //             {
1020 //                 if (run_type == "MC" && evt_NTrack_MC < draw_evt_cut && event_i % print_plot_ratio == 0) {evt_reco_track_gr_ZDCA -> SetPoint(evt_reco_track_gr_ZDCA->GetN(),Get_eta_pair.second, track_phi);}
1021
1022 //                 if (fabs(delta_phi) < phi_diff_cut)
1023 //                 {
1024 //                     if (run_type == "MC" && evt_NTrack_MC < draw_evt_cut && event_i % print_plot_ratio == 0) {evt_reco_track_gr_ZDCAPhi -> SetPoint(evt_reco_track_gr_ZDCAPhi->GetN(),Get_eta_pair.second, track_phi);}
1025 //                     // cout<<"reco eta : "<<Get_eta_pair.second<<" reco phi : "<<track_phi<<endl;
1026
1027 //                     // if (fabs(Get_eta_pair.second) < 0.1) { cout<<"eta "<<Get_eta_pair.second<<", three points : "<<get_radius(beam_origin.first,beam_origin.second)<<" "<<evt_z.first<<" "<<evt_z.second<<" "<<
1028 //                     // 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<<" "<<
1029 //                     // 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<<endl;}
1030
1031 //                     Get_eta_pair = Get_eta(
1032 //                         {get_radius(beam_origin.first,beam_origin.second), evt_z.first,evt_z.second},
1033 //                         {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},
1034 //                         {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}
1035 //                     );
1036
1037 //                     double eta_bin = eta_region_hist -> Fill(Get_eta_pair.second);
1038 //                     // cout<<"test1 "<<eta_bin<<endl;
1039 //                     if (eta_bin != -1)
1040 //                     {
1041 //                         final_track_delta_phi_1D[centrality_map[centrality_bin]][eta_bin - 1]      -> Fill(delta_phi);
1042 //                         final_track_delta_phi_1D[final_track_delta_phi_1D.size() - 1][eta_bin - 1] -> Fill(delta_phi);
1043 //                         out_track_delta_phi_d.push_back(delta_phi);
1044 //                         out_recotrack_eta_d.push_back(Get_eta_pair.second);
1045 //                         out_track_eta_i.push_back(eta_bin);
1046 //                     }
1047
1048
1049 //                     reco_eta_correlation_2D -> Fill(Get_eta_pair.second, (inner_clu_eta + outer_clu_eta)/2.);
1050 //                     reco_eta_diff_reco3P_2D -> Fill(Get_eta_pair.second, Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.);
1051 //                     reco_eta_diff_1D        -> Fill(Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.);
1052
1053 //                     track_DeltaPhi_eta_2D[centrality_map[centrality_bin]]   -> Fill(delta_phi, Get_eta_pair.second);
1054 //                     track_DeltaPhi_eta_2D[track_DeltaPhi_eta_2D.size() - 1] -> Fill(delta_phi, Get_eta_pair.second);
1055
1056 //                     // cout<<"test_5"<<endl;
1057 //                     dNdeta_1D[centrality_map[centrality_bin]] -> Fill(Get_eta_pair.second);
1058 //                     dNdeta_1D[dNdeta_1D.size() - 1]           -> Fill(Get_eta_pair.second);
1059
1060 //                     track_eta_phi_2D[centrality_map[centrality_bin]] -> Fill(track_phi, Get_eta_pair.second);
1061 //                     track_eta_phi_2D[track_eta_phi_2D.size() - 1]    -> Fill(track_phi, Get_eta_pair.second);
1062
1063 //                     track_eta_z_2D[centrality_map[centrality_bin]] -> Fill(Get_eta_pair.second, evt_z.first);
1064 //                     track_eta_z_2D[track_eta_z_2D.size() - 1]      -> Fill(Get_eta_pair.second, evt_z.first);
1065
1066 //                     evt_NTrack += 1;
1067 //                     // cout<<event_i<<" test_10 "<<inner_i<<" "<<outer_i<<endl;
1068 //                 }
1069
1070 //             }
1071
1072 //         }
1073 //     }
1074
1075 //     // note : prepare the tree out
1076 //     out_eID = event_i;
1077 //     out_evt_centrality_bin = centrality_bin;
1078 //     out_evt_zvtx = evt_z.first;
1079 //     tree_out -> Fill();
1080
1081 //     // cout<<" "<<endl;
1082 //     // cout<<" "<<endl;
1083 //     // cout<<"test_8"<<endl;
1084 //     if (run_type == "MC")
1085 //     {
1086 //         track_correlation_2D -> Fill(evt_NTrack_MC, evt_NTrack);
1087 //         track_ratio_1D[centrality_map[centrality_bin]] -> Fill( double(evt_NTrack) / double(evt_NTrack_MC) );
1088 //         track_ratio_1D[track_ratio_1D.size() - 1]      -> Fill( double(evt_NTrack) / double(evt_NTrack_MC) );
1089 //     }
1090
1091 //     track_cluster_ratio_multiplicity_2D -> Fill( total_NClus, double(total_NClus) / double(evt_NTrack) );
1092 //     track_cluster_ratio_1D[centrality_map[centrality_bin]]    -> Fill( double(total_NClus) / double(evt_NTrack) );
1093 //     track_cluster_ratio_1D[track_cluster_ratio_1D.size() - 1] -> Fill( double(total_NClus) / double(evt_NTrack) );
1094
1095 //     if (run_type == "MC" && evt_NTrack_MC < draw_evt_cut && event_i % print_plot_ratio == 0){print_evt_plot(event_i, evt_NTrack_MC, inner_NClu, outer_NClu);}
1096
1097 //     return_tag = 1;
1098 // }
1099
1100 // note : this function is for the event by event vertex calculation
1101 // note : this function is the new method, which means that we first put the cluster into a phi map, and then there are two for loops still, but we only check the certain range of the phi
1102 // note : this function requires that one cluster can only be used once for single tracklet
1103 // note : so maybe there should be no background subtraction for this case
1104 void INTTEta::ProcessEvt(
1105     int event_i,
1106     vector<clu_info> temp_sPH_inner_nocolumn_vec,
1107     vector<clu_info> temp_sPH_outer_nocolumn_vec,
1108     vector<vector<double>> temp_sPH_nocolumn_vec,
1109     vector<vector<double>> temp_sPH_nocolumn_rz_vec,
1110     int NvtxMC,
1111     vector<double> TrigvtxMC,
1112     bool PhiCheckTag,
1113     Long64_t bco_full,
1114     pair<double,double> evt_z,
1115     int centrality_bin,
1116     vector<vector<float>> true_track_info,
1117     double zvtx_weighting = 1.0
1118 )
1119 { // note : evt_z : {z, width}
1120     return_tag = 0;
1121
1122     if (event_i%100 == 0) {cout<<"In INTTEta class, running event : "<<event_i<<endl;}
1123
1124     inner_NClu = temp_sPH_inner_nocolumn_vec.size();
1125     outer_NClu = temp_sPH_outer_nocolumn_vec.size();
1126     total_NClus = inner_NClu + outer_NClu;
1127
1128     // cout<<"test_0"<<endl;
1129     if (total_NClus < zvtx_cal_require) {return; cout<<"return confirmation"<<endl;}
1130
1131     if (run_type == "MC" && NvtxMC != 1) { return; cout<<"In INTTEta class, event : "<<event_i<<" Nvtx : "<<NvtxMC<<" Nvtx more than one "<<endl;}
1132     if (PhiCheckTag == false)            { return; cout<<"In INTTEta class, event : "<<event_i<<" Nvtx : "<<NvtxMC<<" Not full phi has hits "<<endl;}
1133
1134     if (inner_NClu < 10 || outer_NClu < 10 || total_NClus > N_clu_cut || total_NClus < N_clu_cutl)
1135     {
1136         return;
1137         printf("In INTTEta class, event : %i, low clu continue, NClus : %lu \n", event_i, total_NClus);
1138     }
1139
1140     // todo : the z vertex range is here
1141     // if (GetTH2BinXY(coarse_eta_z_map->GetNbinsX(), coarse_eta_z_map->GetNbinsY(), coarse_eta_z_map -> Fill(0., evt_z.first)).second == -1) {return;}
1142     // if (-220 > evt_z.first + evt_z.second || -180 < evt_z.first - evt_z.second) {return;}
1143     // if (-100 > evt_z.first + evt_z.second || 100 < evt_z.first - evt_z.second) {return;}
1144
1145     // note : for data, the denominator
1146     MBin_Z_evt_hist_2D -> Fill(centrality_map[centrality_bin], evt_z.first, zvtx_weighting);
1147     MBin_Z_evt_hist_2D -> Fill(MBin_Z_evt_hist_2D -> GetNbinsX()-1, evt_z.first, zvtx_weighting);
1148
1149     N_GoodEvent += 1;
1150     N_GoodEvent_vec[centrality_map[centrality_bin]] += 1;
1151
1152     // cout<<"N inner cluster : "<<inner_NClu<<" N outer cluster : "<<outer_NClu<<endl;
1153
1154     double INTT_eta_acceptance_l = -0.5 * TMath::Log((sqrt(pow(-230.-TrigvtxMC[2]*10.,2)+pow(INTT_layer_R[3],2))-(-230.-TrigvtxMC[2]*10.)) / (sqrt(pow(-230.-TrigvtxMC[2]*10.,2)+pow(INTT_layer_R[3],2))+(-230.-TrigvtxMC[2]*10.))); // note : left
1155     double INTT_eta_acceptance_r =  -0.5 * TMath::Log((sqrt(pow(230.-TrigvtxMC[2]*10.,2)+pow(INTT_layer_R[3],2))-(230.-TrigvtxMC[2]*10.)) / (sqrt(pow(230.-TrigvtxMC[2]*10.,2)+pow(INTT_layer_R[3],2))+(230.-TrigvtxMC[2]*10.))); // note : right
1156
1157     if (run_type == "MC")
1158     {
1159         MBin_Z_evt_hist_2D_MC -> Fill(centrality_map[centrality_bin], TrigvtxMC[2]*10.);
1160         MBin_Z_evt_hist_2D_MC -> Fill(MBin_Z_evt_hist_2D_MC -> GetNbinsX()-1, TrigvtxMC[2]*10.);
1161
1162         // note : for the true track case ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1163         // if (event_i % 100 == 0){cout<<"z : "<<TrigvtxMC[2]*10.<<" eta : "<<INTT_eta_acceptance_l<<" "<<INTT_eta_acceptance_r<<endl;}
1164
1165         // cout<<"true_track_info : "<<true_track_info.size()<<endl;
1166         for (int track_i = 0; track_i < true_track_info.size(); track_i++)
1167         {
1168             if (true_track_info[track_i][2] == 111 || true_track_info[track_i][2] == 22 || abs(true_track_info[track_i][2]) == 2112){continue;}
1169
1170             // note : coarse eta and z binning, this is for the final dNdeta result
1171             coarse_eta_z_2D_fulleta_MC[centrality_map[centrality_bin]]        -> Fill(true_track_info[track_i][0], TrigvtxMC[2]*10.);
1172             coarse_eta_z_2D_fulleta_MC[coarse_eta_z_2D_fulleta_MC.size() - 1] -> Fill(true_track_info[track_i][0], TrigvtxMC[2]*10.);
1173
1174             if (TrigvtxMC[2]*10. > z_region[selected_z_region_id.first] && TrigvtxMC[2]*10. < z_region[selected_z_region_id.second])
1175             {
1176                 dNdeta_1D_MC[centrality_map[centrality_bin]] -> Fill(true_track_info[track_i][0]);
1177                 dNdeta_1D_MC[dNdeta_1D_MC.size() - 1]        -> Fill(true_track_info[track_i][0]);
1178             }
1179
1180
1181             if (true_track_info[track_i][0] > INTT_eta_acceptance_l && true_track_info[track_i][0] < INTT_eta_acceptance_r)
1182             {
1183                 // note : 1D, fine binning for a detailed check, full zvtx range
1184                 if (TrigvtxMC[2]*10. > z_region[selected_z_region_id.first] && TrigvtxMC[2]*10. < z_region[selected_z_region_id.second])
1185                 {
1186                     dNdeta_1D_MC_edge_eta_cut[centrality_map[centrality_bin]]       -> Fill(true_track_info[track_i][0]);
1187                     dNdeta_1D_MC_edge_eta_cut[dNdeta_1D_MC_edge_eta_cut.size() - 1] -> Fill(true_track_info[track_i][0]);
1188                 }
1189
1190                 if (GetTH2BinXY(coarse_eta_z_map -> GetNbinsX(), coarse_eta_z_map -> GetNbinsY(), coarse_eta_z_map -> Fill(true_track_info[track_i][0], TrigvtxMC[2]*10.)).second == tight_zvtx_bin)
1191                 {
1192                     // note : 1D, coarse binning for the final dNdeta result with a fixed z vertex range
1193                     final_dNdeta_1D_MC[centrality_map[centrality_bin]] -> Fill(true_track_info[track_i][0]);
1194                     final_dNdeta_1D_MC[final_dNdeta_1D_MC.size() - 1]  -> Fill(true_track_info[track_i][0]);
1195                 }
1196
1197                 // note : 2D, fine binning for the detailed check of eta and z distribution
1198                 track_eta_z_2D_MC[centrality_map[centrality_bin]] -> Fill(true_track_info[track_i][0], TrigvtxMC[2]*10.);
1199                 track_eta_z_2D_MC[track_eta_z_2D_MC.size() - 1]   -> Fill(true_track_info[track_i][0], TrigvtxMC[2]*10.);
1200
1201                 // note : coarse eta and z binning, this is for the final dNdeta result
1202                 coarse_eta_z_2D_MC[centrality_map[centrality_bin]] -> Fill(true_track_info[track_i][0], TrigvtxMC[2]*10.);
1203                 coarse_eta_z_2D_MC[coarse_eta_z_2D_MC.size() - 1]  -> Fill(true_track_info[track_i][0], TrigvtxMC[2]*10.);
1204
1205
1206                 // cout<<"true track eta : "<<true_track_info[track_i][0]<<" phi : "<<convertTo360(true_track_info[track_i][1])<<endl;
1207                 // cout<<"("<<true_track_info[track_i][0]<<", "<<convertTo360(true_track_info[track_i][1])<<"), ";
1208
1209                 evt_true_track_gr -> SetPoint(evt_true_track_gr->GetN(),true_track_info[track_i][0], convertTo360(true_track_info[track_i][1]));
1210                 out_truetrack_eta_d.push_back(true_track_info[track_i][0]);
1211                 out_truetrack_phi_d.push_back(convertTo360(true_track_info[track_i][1]));
1212
1213                 evt_NTrack_MC += 1;
1214             }
1215         }
1216         // if (evt_NTrack_MC < 10) {cout<<"evt : "<<event_i<<" ---- N reco track : "<<evt_NTrack<<" N true track : "<<evt_NTrack_MC<<" ratio : "<<double(evt_NTrack) / double(evt_NTrack_MC)<<endl;}
1217     }
1218
1219     if (evt_z.first > z_region[selected_z_region_id.first] && evt_z.first < z_region[selected_z_region_id.second])
1220     {
1221         exclusive_NClus_inner -> Fill(temp_sPH_inner_nocolumn_vec.size(), zvtx_weighting);
1222         exclusive_NClus_outer -> Fill(temp_sPH_outer_nocolumn_vec.size(), zvtx_weighting);
1223         exclusive_NClus_sum   -> Fill(temp_sPH_inner_nocolumn_vec.size() + temp_sPH_outer_nocolumn_vec.size(), zvtx_weighting);
1224     }
1225
1226     // note : put the cluster into the phi map, the first bool is for the cluster usage.
1227     // note : false means the cluster is not used
1228     for (int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++) {
1229         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());
1230         double Clus_InnerPhi_Offset_radian = atan2(temp_sPH_inner_nocolumn_vec[inner_i].y - beam_origin.second, temp_sPH_inner_nocolumn_vec[inner_i].x - beam_origin.first);
1231
1232         // cout<<"inner clu phi : "<<Clus_InnerPhi_Offset<<" origin: "<< temp_sPH_inner_nocolumn_vec[inner_i].phi <<endl;
1233         // cout<<" ("<<Clus_InnerPhi_Offset<<", "<< temp_sPH_inner_nocolumn_vec[inner_i].phi<<")" <<endl;
1234         inner_clu_phi_map[ int(Clus_InnerPhi_Offset) ].push_back({false,temp_sPH_inner_nocolumn_vec[inner_i]});
1235         check_inner_layer_clu_phi_1D->Fill(Clus_InnerPhi_Offset);
1236
1237         double clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y, temp_sPH_inner_nocolumn_vec[inner_i].z});
1238
1239         if (evt_z.first > z_region[selected_z_region_id.first] && evt_z.first < z_region[selected_z_region_id.second]){
1240             exclusive_cluster_inner_eta -> Fill(clu_eta, zvtx_weighting);
1241             exclusive_cluster_inner_phi -> Fill(Clus_InnerPhi_Offset_radian, zvtx_weighting);
1242
1243             exclusive_cluster_inner_eta_phi_2D -> Fill(clu_eta, Clus_InnerPhi_Offset_radian, zvtx_weighting);
1244
1245             exclusive_cluster_inner_eta_adc_2D -> Fill(clu_eta, temp_sPH_inner_nocolumn_vec[inner_i].sum_adc, zvtx_weighting);
1246
1247             exclusive_cluster_all_eta   -> Fill(clu_eta, zvtx_weighting);
1248             exclusive_cluster_all_phi   -> Fill(Clus_InnerPhi_Offset_radian, zvtx_weighting);
1249
1250             exclusive_cluster_inner_adc -> Fill(temp_sPH_inner_nocolumn_vec[inner_i].sum_adc, zvtx_weighting);
1251         }
1252
1253         if (clu_eta > INTT_eta_acceptance_l && clu_eta < INTT_eta_acceptance_r) {effective_total_NClus += 1;}
1254     }
1255     for (int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++) {
1256         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());
1257         double Clus_OuterPhi_Offset_radian = atan2(temp_sPH_outer_nocolumn_vec[outer_i].y - beam_origin.second, temp_sPH_outer_nocolumn_vec[outer_i].x - beam_origin.first);
1258
1259         outer_clu_phi_map[ int(Clus_OuterPhi_Offset) ].push_back({false,temp_sPH_outer_nocolumn_vec[outer_i]});
1260         check_outer_layer_clu_phi_1D->Fill(Clus_OuterPhi_Offset);
1261
1262         double clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y, temp_sPH_outer_nocolumn_vec[outer_i].z});
1263
1264         if (evt_z.first > z_region[selected_z_region_id.first] && evt_z.first < z_region[selected_z_region_id.second]){
1265             exclusive_cluster_outer_eta -> Fill(clu_eta, zvtx_weighting);
1266             exclusive_cluster_outer_phi -> Fill(Clus_OuterPhi_Offset_radian, zvtx_weighting);
1267
1268             exclusive_cluster_outer_eta_phi_2D -> Fill(clu_eta, Clus_OuterPhi_Offset_radian, zvtx_weighting);
1269
1270             exclusive_cluster_outer_eta_adc_2D -> Fill(clu_eta, temp_sPH_outer_nocolumn_vec[outer_i].sum_adc, zvtx_weighting);
1271
1272             exclusive_cluster_all_eta   -> Fill(clu_eta, zvtx_weighting);
1273             exclusive_cluster_all_phi   -> Fill(Clus_OuterPhi_Offset_radian, zvtx_weighting);
1274
1275             exclusive_cluster_outer_adc -> Fill(temp_sPH_outer_nocolumn_vec[outer_i].sum_adc, zvtx_weighting);
1276         }
1277
1278         if (clu_eta > INTT_eta_acceptance_l && clu_eta < INTT_eta_acceptance_r) {effective_total_NClus += 1;}
1279     }
1280
1281     // note : for the mega cluster finder, the 3/4-cluster tracklet that happens at the overlap region
1282     // mega_track_finder -> FindMegaTracks(inner_clu_phi_map, outer_clu_phi_map, evt_z, centrality_map[centrality_bin]);
1283     // note : the followings are the study of the mega track finder efficiency
1284     // int N_mega_track_with_recoZ = mega_track_finder -> Get_NClu3_track_count() + mega_track_finder -> Get_NClu4_track_count();
1285     // mega_track_finder -> ClearEvt();
1286     // mega_track_finder -> FindMegaTracks(inner_clu_phi_map, outer_clu_phi_map, {TrigvtxMC[2] * 10., 0.5}, centrality_map[centrality_bin]);
1287     // int N_mega_track_with_trueZ = mega_track_finder -> Get_NClu3_track_count() + mega_track_finder -> Get_NClu4_track_count();
1288     // mega_track_finder -> ClearEvt();
1289     // if (N_mega_track_with_recoZ != N_mega_track_with_trueZ) {cout<<N_mega_track_with_recoZ<<" "<<N_mega_track_with_trueZ<<endl;}
1290     // if (N_mega_track_with_trueZ != 0) {mega_track_finding_ratio_2D -> Fill(total_NClus, double(N_mega_track_with_recoZ) / double(N_mega_track_with_trueZ));}
1291
1292     // note : for two-cluster tracklets only
1293     for (int inner_phi_i = 0; inner_phi_i < 360; inner_phi_i++) // note : each phi cell (1 degree)
1294     {
1295         // note : N cluster in this phi cell
1296         for (int inner_phi_clu_i = 0; inner_phi_clu_i < inner_clu_phi_map[inner_phi_i].size(); inner_phi_clu_i++)
1297         {
1298             clu_multi_used_tight = 0;
1299             clu_multi_used_loose = 0;
1300             if (inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].first == true) {continue;}
1301
1302             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());
1303             double Clus_InnerPhi_Offset_radian = 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);
1304
1305             // todo: change the outer phi scan range
1306             // note : the outer phi index, -4, -3, -2, -1, 0, 1, 2, 3, 4
1307             for (int scan_i = -4; scan_i < 5; scan_i++)
1308             {
1309                 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;
1310
1311                 // note : N clusters in that outer phi cell
1312                 for (int outer_phi_clu_i = 0; outer_phi_clu_i < outer_clu_phi_map[true_scan_i].size(); outer_phi_clu_i++)
1313                 {
1314                     if (outer_clu_phi_map[true_scan_i][outer_phi_clu_i].first == true) {continue;}
1315
1316                     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());
1317                     double Clus_OuterPhi_Offset_radian = 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);
1318                     double delta_phi = get_delta_phi(Clus_InnerPhi_Offset, Clus_OuterPhi_Offset);
1319
1320                     // if (fabs(delta_phi) > 5.72) {continue;}
1321                     if (fabs(delta_phi) > 3.5) {continue;}
1322
1323                     double inner_clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{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, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z});
1324                     double outer_clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{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, outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.z});
1325                     double delta_eta = inner_clu_eta - outer_clu_eta;
1326
1327                     // note : all zvtx range and all eta region
1328                     track_delta_eta_1D[centrality_map[centrality_bin]] -> Fill( delta_eta );
1329                     track_delta_eta_1D[track_delta_eta_1D.size() - 1]  -> Fill( delta_eta );
1330
1331                     // note : all zvtx range and all eta region
1332                     track_DeltaPhi_DeltaEta_2D[centrality_map[centrality_bin]]        -> Fill(delta_phi, delta_eta);
1333                     track_DeltaPhi_DeltaEta_2D[track_DeltaPhi_DeltaEta_2D.size() - 1] -> Fill(delta_phi, delta_eta);
1334
1335                     pair<double,double> z_range_info = Get_possible_zvtx(
1336                         0., // get_radius(beam_origin.first,beam_origin.second),
1337                         {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
1338                         {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
1339                     );
1340
1341                     // 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
1342                     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;}
1343                     if (fabs(delta_phi) < 0.6) {clu_multi_used_tight += 1;} // todo : the delta phi cut value of the multi-usage inner cluster is here
1344                     track_delta_eta_1D_post[centrality_map[centrality_bin]]     -> Fill( delta_eta );
1345                     track_delta_eta_1D_post[track_delta_eta_1D_post.size() - 1] -> Fill( delta_eta );
1346
1347                     double DCA_sign = calculateAngleBetweenVectors(
1348                         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,
1349                         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,
1350                         beam_origin.first, beam_origin.second
1351                     );
1352
1353                     track_DCA_distance[centrality_map[centrality_bin]] -> Fill( DCA_sign );
1354                     track_DCA_distance[track_DCA_distance.size() - 1]  -> Fill( DCA_sign );
1355
1356                     track_phi_DCA_2D[centrality_map[centrality_bin]] -> Fill( delta_phi, DCA_sign );
1357                     track_phi_DCA_2D[track_phi_DCA_2D.size() - 1]    -> Fill( delta_phi, DCA_sign );
1358
1359                     track_delta_phi_1D[centrality_map[centrality_bin]] -> Fill( delta_phi ); // note : the pairs that pass the rough delta phi cut and eta cut
1360                     track_delta_phi_1D[track_delta_phi_1D.size() - 1]  -> Fill( delta_phi );
1361
1362                     proto_pair_index.push_back({inner_phi_i, inner_phi_clu_i, true_scan_i, outer_phi_clu_i});
1363                     proto_pair_delta_phi_abs.push_back(fabs(delta_phi));
1364                     proto_pair_delta_phi.push_back({Clus_InnerPhi_Offset, Clus_OuterPhi_Offset, delta_phi});
1365
1366
1367                     // note : do the fill here (find the best match outer cluster with the inner cluster )
1368                     Get_eta_pair = Get_eta(
1369                         {0., evt_z.first,evt_z.second},
1370                         {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},
1371                         {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}
1372                     );
1373
1374                     // note : find the bin of the eta
1375                     // double eta_bin = eta_region_hist -> Fill(Get_eta_pair.second);
1376                     // // cout<<"test1 "<<eta_bin<<endl;
1377                     // if (eta_bin != -1)
1378                     // {
1379                     //     final_track_multi_delta_phi_1D[centrality_map[centrality_bin]][eta_bin - 1]            -> Fill(delta_phi);
1380                     //     final_track_multi_delta_phi_1D[final_track_multi_delta_phi_1D.size() - 1][eta_bin - 1] -> Fill(delta_phi);
1381
1382                     //     if (fabs(delta_phi) <= signal_region) {
1383                     //         good_tracklet_multi_counting[centrality_map[centrality_bin]][eta_bin - 1]          += 1;
1384                     //         good_tracklet_multi_counting[good_tracklet_multi_counting.size() - 1][eta_bin - 1] += 1;
1385                     //     }
1386                     // }
1387
1388                     // note : the code below is for the 2D eta-z array
1389                     double eta_z_bin = GetTH2Index1D(GetTH2BinXY(coarse_eta_z_map->GetNbinsX(), coarse_eta_z_map->GetNbinsY(), coarse_eta_z_map -> Fill(Get_eta_pair.second, evt_z.first)), coarse_eta_z_map->GetNbinsX());
1390                     if (int(eta_z_bin) != -1)
1391                     {
1392                         final_track_multi_delta_phi_1D[centrality_map[centrality_bin]][eta_z_bin]            -> Fill(delta_phi, zvtx_weighting); // note : each coarse eta and z bin
1393                         final_track_multi_delta_phi_1D[final_track_multi_delta_phi_1D.size() - 1][eta_z_bin] -> Fill(delta_phi, zvtx_weighting);
1394
1395                         if (fabs(delta_phi) <= signal_region) {
1396                             good_tracklet_multi_counting[centrality_map[centrality_bin]][eta_z_bin]          += 1;
1397                             good_tracklet_multi_counting[good_tracklet_multi_counting.size() - 1][eta_z_bin] += 1;
1398
1399                             coarse_Reco_SignalNTracklet_Multi_eta_z_2D[centrality_map[centrality_bin]]                        -> Fill(Get_eta_pair.second, evt_z.first, zvtx_weighting);
1400                             coarse_Reco_SignalNTracklet_Multi_eta_z_2D[coarse_Reco_SignalNTracklet_Multi_eta_z_2D.size() - 1] -> Fill(Get_eta_pair.second, evt_z.first, zvtx_weighting);
1401
1402                             if (evt_z.first > z_region[selected_z_region_id.first] && evt_z.first < z_region[selected_z_region_id.second])
1403                             {
1404                                 exclusive_loose_tracklet_eta -> Fill(Get_eta_pair.second, zvtx_weighting);
1405                                 exclusive_loose_tracklet_phi -> Fill( (Clus_InnerPhi_Offset_radian + Clus_OuterPhi_Offset_radian)/2. , zvtx_weighting);
1406                             }
1407
1408                         }
1409
1410                     }
1411
1412                     // if (clu_multi_used_loose == 0)
1413                     // {
1414                     //     pair_delta_phi = delta_phi;
1415                     //     pair_outer_index = {true_scan_i, outer_phi_clu_i};
1416                     //     clu_multi_used_loose += 1;
1417                     // }
1418                     // else if (fabs(delta_phi) < fabs(pair_delta_phi))
1419                     // {
1420                     //     pair_delta_phi = delta_phi;
1421                     //     pair_outer_index = {true_scan_i, outer_phi_clu_i};
1422                     //     clu_multi_used_loose += 1;
1423                     // }
1424                 }
1425             } // note : end outer loop
1426
1427             clu_used_centrality_2D -> Fill(total_NClus, clu_multi_used_tight);
1428
1429             // // note : if there is no good match with inner and outer clusters, continue
1430             // if (clu_multi_used_loose == 0) {continue;}
1431
1432             // // note : do the fill here (find the best match outer cluster with the inner cluster )
1433             // Get_eta_pair = Get_eta(
1434             //     {0., evt_z.first,evt_z.second},
1435             //     {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},
1436             //     {get_radius(outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.x - beam_origin.first, outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.y - beam_origin.second), outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.z}
1437             // );
1438
1439             // double inner_clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{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, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z});
1440             // double outer_clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.x, outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.y, outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.z});
1441             // double track_phi = Clus_InnerPhi_Offset - (pair_delta_phi/2.);
1442
1443             // if  (Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2. > 0.3)
1444             // {
1445             //     cout<<" "<<endl;
1446             //     cout<<"inner clu eta : "<<inner_clu_eta<<" outer clu eta : "<<outer_clu_eta<<"avg eta : "<< (inner_clu_eta + outer_clu_eta)/2. <<" reco eta : "<<Get_eta_pair.second<<" diff: "<<Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.<<endl;
1447             //     cout<<"inner clu pos : "<<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<<" "<<inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z<<endl;
1448             //     cout<<"outer clu pos : "<<outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.x<<" "<<outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.y<<" "<<outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.z<<endl;
1449             // }
1450
1451             // reco_eta_correlation_2D -> Fill(Get_eta_pair.second, (inner_clu_eta + outer_clu_eta)/2.);
1452             // reco_eta_diff_reco3P_2D -> Fill(Get_eta_pair.second, Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.);
1453             // reco_eta_diff_1D        -> Fill(Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.);
1454
1455             // // track_delta_phi_1D[centrality_map[centrality_bin]] -> Fill( pair_delta_phi );
1456             // // track_delta_phi_1D[track_delta_phi_1D.size() - 1]  -> Fill( pair_delta_phi );
1457
1458             // track_DeltaPhi_eta_2D[centrality_map[centrality_bin]]   -> Fill(pair_delta_phi, Get_eta_pair.second);
1459             // track_DeltaPhi_eta_2D[track_DeltaPhi_eta_2D.size() - 1] -> Fill(pair_delta_phi, Get_eta_pair.second);
1460
1461             // // cout<<"test_5"<<endl;
1462             // dNdeta_1D[centrality_map[centrality_bin]] -> Fill(Get_eta_pair.second);
1463             // dNdeta_1D[dNdeta_1D.size() - 1]           -> Fill(Get_eta_pair.second);
1464
1465             // track_eta_phi_2D[centrality_map[centrality_bin]] -> Fill(track_phi, Get_eta_pair.second);
1466             // track_eta_phi_2D[track_eta_phi_2D.size() - 1]    -> Fill(track_phi, Get_eta_pair.second);
1467
1468             // track_eta_z_2D[centrality_map[centrality_bin]] -> Fill(Get_eta_pair.second, evt_z.first);
1469             // track_eta_z_2D[track_eta_z_2D.size() - 1]      -> Fill(Get_eta_pair.second, evt_z.first);
1470
1471             // // note : find the bin of the eta
1472             // double eta_bin = eta_region_hist -> Fill(Get_eta_pair.second);
1473             // // cout<<"test1 "<<eta_bin<<endl;
1474             // if (eta_bin != -1)
1475             // {
1476             //     final_track_delta_phi_1D[centrality_map[centrality_bin]][eta_bin - 1]      -> Fill(pair_delta_phi);
1477             //     final_track_delta_phi_1D[final_track_delta_phi_1D.size() - 1][eta_bin - 1] -> Fill(pair_delta_phi);
1478
1479             //     out_track_delta_phi_d.push_back(pair_delta_phi);
1480             //     out_recotrack_eta_d.push_back(Get_eta_pair.second);
1481             //     out_track_eta_i.push_back(eta_bin);
1482             // }
1483
1484             // evt_NTrack += 1;
1485             // // note : this outer cluster is used
1486             // outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].first = true;
1487         }
1488     } // note : end inner loop
1489
1490     // note : if there is no good proto-track then move on to the next event
1491     if (proto_pair_delta_phi_abs.size() == 0) {return;}
1492
1493     long long vec_size = proto_pair_delta_phi_abs.size();
1494     long long ind[proto_pair_delta_phi_abs.size()];
1495     TMath::Sort(vec_size, &proto_pair_delta_phi_abs[0], ind, false);
1496     for (int pair_i = 0; pair_i < proto_pair_index.size(); pair_i++)
1497     {
1498         int inner_index_0 = proto_pair_index[ind[pair_i]][0];
1499         int inner_index_1 = proto_pair_index[ind[pair_i]][1];
1500         int outer_index_0 = proto_pair_index[ind[pair_i]][2];
1501         int outer_index_1 = proto_pair_index[ind[pair_i]][3];
1502
1503         // cout<<"test: the pair delta phi abs: "<<proto_pair_delta_phi_abs[ind[pair_i]]<<endl;
1504         if (inner_used_clu[Form("%i_%i", inner_index_0, inner_index_1)] != 0) { inner_used_clu[Form("%i_%i", inner_index_0, inner_index_1)] += 1; continue;}
1505         if (outer_used_clu[Form("%i_%i", outer_index_0, outer_index_1)] != 0) { outer_used_clu[Form("%i_%i", outer_index_0, outer_index_1)] += 1; continue;}
1506         double delta_phi = proto_pair_delta_phi[ind[pair_i]][2];
1507         double track_phi = get_track_phi(proto_pair_delta_phi[ind[pair_i]][0], delta_phi); // proto_pair_delta_phi[ind[pair_i]][0] - (delta_phi/2.);
1508
1509         // note : do the fill here (find the best match outer cluster with the inner cluster )
1510         Get_eta_pair = Get_eta(
1511             {0., evt_z.first,evt_z.second},
1512             {get_radius(inner_clu_phi_map[inner_index_0][inner_index_1].second.x - beam_origin.first, inner_clu_phi_map[inner_index_0][inner_index_1].second.y - beam_origin.second), inner_clu_phi_map[inner_index_0][inner_index_1].second.z},
1513             {get_radius(outer_clu_phi_map[outer_index_0][outer_index_1].second.x - beam_origin.first, outer_clu_phi_map[outer_index_0][outer_index_1].second.y - beam_origin.second), outer_clu_phi_map[outer_index_0][outer_index_1].second.z}
1514         );
1515
1516         double inner_clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{inner_clu_phi_map[inner_index_0][inner_index_1].second.x, inner_clu_phi_map[inner_index_0][inner_index_1].second.y, inner_clu_phi_map[inner_index_0][inner_index_1].second.z});
1517         double outer_clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{outer_clu_phi_map[outer_index_0][outer_index_1].second.x, outer_clu_phi_map[outer_index_0][outer_index_1].second.y, outer_clu_phi_map[outer_index_0][outer_index_1].second.z});
1518
1519         double Clus_InnerPhi_Offset_radian = atan2(inner_clu_phi_map[inner_index_0][inner_index_1].second.y - beam_origin.second, inner_clu_phi_map[inner_index_0][inner_index_1].second.x - beam_origin.first);
1520         double Clus_OuterPhi_Offset_radian = atan2(outer_clu_phi_map[outer_index_0][outer_index_1].second.y - beam_origin.second, outer_clu_phi_map[outer_index_0][outer_index_1].second.x - beam_origin.first);
1521
1522         if  (Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2. > 0.3)
1523         {
1524             cout<<" "<<endl;
1525             cout<<"inner clu eta : "<<inner_clu_eta<<" outer clu eta : "<<outer_clu_eta<<"avg eta : "<< (inner_clu_eta + outer_clu_eta)/2. <<" reco eta : "<<Get_eta_pair.second<<" diff: "<<Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.<<endl;
1526             cout<<"inner clu pos : "<<inner_clu_phi_map[inner_index_0][inner_index_1].second.x<<" "<<inner_clu_phi_map[inner_index_0][inner_index_1].second.y<<" "<<inner_clu_phi_map[inner_index_0][inner_index_1].second.z<<endl;
1527             cout<<"outer clu pos : "<<outer_clu_phi_map[outer_index_0][outer_index_1].second.x<<" "<<outer_clu_phi_map[outer_index_0][outer_index_1].second.y<<" "<<outer_clu_phi_map[outer_index_0][outer_index_1].second.z<<endl;
1528         }
1529
1530         reco_eta_correlation_2D -> Fill(Get_eta_pair.second, (inner_clu_eta + outer_clu_eta)/2.);
1531         reco_eta_diff_reco3P_2D -> Fill(Get_eta_pair.second, Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.);
1532         reco_eta_diff_1D        -> Fill(Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.);
1533
1534         // track_delta_phi_1D[centrality_map[centrality_bin]] -> Fill( delta_phi );
1535         // track_delta_phi_1D[track_delta_phi_1D.size() - 1]  -> Fill( delta_phi );
1536
1537         track_DeltaPhi_eta_2D[centrality_map[centrality_bin]]   -> Fill(delta_phi, Get_eta_pair.second);
1538         track_DeltaPhi_eta_2D[track_DeltaPhi_eta_2D.size() - 1] -> Fill(delta_phi, Get_eta_pair.second);
1539
1540         // cout<<"test_5"<<endl;
1541         if (fabs(delta_phi) <= signal_region)
1542         {
1543             if (evt_z.first > z_region[selected_z_region_id.first] && evt_z.first < z_region[selected_z_region_id.second])
1544             {
1545                 dNdeta_1D[centrality_map[centrality_bin]] -> Fill(Get_eta_pair.second);
1546                 dNdeta_1D[dNdeta_1D.size() - 1]           -> Fill(Get_eta_pair.second);
1547             }
1548
1549
1550             track_eta_z_2D[centrality_map[centrality_bin]] -> Fill(Get_eta_pair.second, evt_z.first);
1551             track_eta_z_2D[track_eta_z_2D.size() - 1]      -> Fill(Get_eta_pair.second, evt_z.first);
1552
1553             out_recotrack_eta_d.push_back(Get_eta_pair.second);
1554             out_recotrack_phi_d.push_back(track_phi);
1555
1556             // cout<<"good 2clu track, eta: "<<Get_eta_pair.second<<" phi: "<<track_phi<<" delta_phi: "<<delta_phi<<" -- index, inner: "<<inner_index_0<<" "<<inner_index_1<<" outer: "<<outer_index_0<<" "<<outer_index_1<<endl;
1557             // cout<<" pair_line->DrawLine("<<inner_index_0<<","<<inner_index_1<<", "<<outer_index_0<<","<<outer_index_1<<");"<<endl;
1558             // cout<<"("<<track_phi<<", "<<Get_eta_pair.second<<"), ";
1559             evt_NTrack += 1;
1560         }
1561
1562
1563         track_eta_phi_2D[centrality_map[centrality_bin]] -> Fill(track_phi, Get_eta_pair.second);
1564         track_eta_phi_2D[track_eta_phi_2D.size() - 1]    -> Fill(track_phi, Get_eta_pair.second);
1565
1566         // // note : find the bin of the eta
1567         // double eta_bin = eta_region_hist -> Fill(Get_eta_pair.second);
1568         // // cout<<"test1 "<<eta_bin<<endl;
1569         // if (eta_bin != -1)
1570         // {
1571         //     final_track_delta_phi_1D[centrality_map[centrality_bin]][eta_bin - 1]      -> Fill(delta_phi);
1572         //     final_track_delta_phi_1D[final_track_delta_phi_1D.size() - 1][eta_bin - 1] -> Fill(delta_phi);
1573
1574         //     if (fabs(delta_phi) <= signal_region) {
1575         //         good_tracklet_counting[centrality_map[centrality_bin]][eta_bin - 1] += 1;
1576         //         good_tracklet_counting[good_tracklet_counting.size() - 1][eta_bin - 1] += 1;
1577         //     }
1578
1579         //     out_track_delta_phi_d.push_back(delta_phi);
1580         //     out_recotrack_eta_d.push_back(Get_eta_pair.second);
1581         //     out_track_eta_i.push_back(eta_bin);
1582         // }
1583
1584         // note : the code below is for the 2D eta-z array
1585         double eta_z_bin = GetTH2Index1D(GetTH2BinXY(coarse_eta_z_map->GetNbinsX(), coarse_eta_z_map->GetNbinsY(), coarse_eta_z_map -> Fill(Get_eta_pair.second, evt_z.first)), coarse_eta_z_map->GetNbinsX());
1586         if (int(eta_z_bin) != -1)
1587         {
1588             final_track_delta_phi_1D[centrality_map[centrality_bin]][eta_z_bin]      -> Fill(delta_phi);
1589             final_track_delta_phi_1D[final_track_delta_phi_1D.size() - 1][eta_z_bin] -> Fill(delta_phi);
1590
1591             if (fabs(delta_phi) <= signal_region) {
1592                 good_tracklet_counting[centrality_map[centrality_bin]][eta_z_bin]    += 1;
1593                 good_tracklet_counting[good_tracklet_counting.size() - 1][eta_z_bin] += 1;
1594
1595                 coarse_Reco_SignalNTracklet_Single_eta_z_2D[centrality_map[centrality_bin]]                         -> Fill(Get_eta_pair.second, evt_z.first, zvtx_weighting);
1596                 coarse_Reco_SignalNTracklet_Single_eta_z_2D[coarse_Reco_SignalNTracklet_Single_eta_z_2D.size() - 1] -> Fill(Get_eta_pair.second, evt_z.first, zvtx_weighting);
1597
1598                 if (evt_z.first > z_region[selected_z_region_id.first] && evt_z.first < z_region[selected_z_region_id.second])
1599                 {
1600                     exclusive_tight_tracklet_eta -> Fill(Get_eta_pair.second, zvtx_weighting);
1601                     exclusive_tight_tracklet_phi -> Fill((Clus_InnerPhi_Offset_radian + Clus_OuterPhi_Offset_radian)/2., zvtx_weighting);
1602                 }
1603             }
1604
1605             // note : save the information in detail
1606             // out_track_delta_phi_d.push_back(delta_phi);
1607             // out_recotrack_eta_d.push_back(Get_eta_pair.second);
1608             // out_track_eta_i.push_back(eta_bin);
1609
1610         }
1611
1612
1613
1614         // note : since the clusters are used in the tracklet, mark the clusters as used
1615         inner_used_clu[Form("%i_%i", inner_index_0, inner_index_1)] += 1;
1616         outer_used_clu[Form("%i_%i", outer_index_0, outer_index_1)] += 1;
1617     }
1618     // note : save the information in detail
1619     out_eID = event_i;
1620     out_evt_centrality_bin = centrality_bin;
1621     out_evt_zvtx = evt_z.first;
1622     out_true_zvtx = TrigvtxMC[2] * 10.;
1623     out_NTrueTrack = evt_NTrack_MC;
1624     out_total_NClus = total_NClus;
1625     out_N2Clu_track = evt_NTrack;
1626     out_N3Clu_track = mega_track_finder -> Get_NClu3_track_count();
1627     out_N4Clu_track = mega_track_finder -> Get_NClu4_track_count();
1628     tree_out -> Fill();
1629
1630     cout<<"event confirmation : "<<out_N2Clu_track<<" "<<out_N3Clu_track<<" "<<out_N4Clu_track<<endl;
1631
1632     // cout<<" "<<endl;
1633     // cout<<" "<<endl;
1634     // cout<<"test_8"<<endl;
1635     if (run_type == "MC")
1636     {
1637         N_reco_cluster_short += (evt_NTrack_MC - evt_NTrack);
1638         track_correlation_2D -> Fill(evt_NTrack_MC, evt_NTrack);
1639         track_ratio_2D -> Fill(evt_NTrack_MC, double(evt_NTrack)/double(evt_NTrack_MC));
1640         track_ratio_1D[centrality_map[centrality_bin]] -> Fill( double(evt_NTrack) / double(evt_NTrack_MC) );
1641         track_ratio_1D[track_ratio_1D.size() - 1]      -> Fill( double(evt_NTrack) / double(evt_NTrack_MC) );
1642
1643         track_cluster_ratio_1D_MC[centrality_map[centrality_bin]]       -> Fill( double(effective_total_NClus) / double(evt_NTrack_MC) );
1644         track_cluster_ratio_1D_MC[track_cluster_ratio_1D_MC.size() - 1] -> Fill( double(effective_total_NClus) / double(evt_NTrack_MC) );
1645         track_cluster_ratio_multiplicity_2D_MC -> Fill( effective_total_NClus, double(effective_total_NClus) / double(evt_NTrack_MC) );
1646     }
1647
1648     track_cluster_ratio_multiplicity_2D -> Fill( effective_total_NClus, double(effective_total_NClus) / double(evt_NTrack) );
1649     track_cluster_ratio_1D[centrality_map[centrality_bin]]    -> Fill( double(effective_total_NClus) / double(evt_NTrack) );
1650     track_cluster_ratio_1D[track_cluster_ratio_1D.size() - 1] -> Fill( double(effective_total_NClus) / double(evt_NTrack) );
1651
1652     // if (run_type == "MC" && evt_NTrack_MC < draw_evt_cut && event_i % print_plot_ratio == 0){print_evt_plot(event_i, evt_NTrack_MC, inner_NClu, outer_NClu);}
1653     return_tag = 1;
1654 }
1655
1656
1657 void INTTEta::ClearEvt()
1658 {
1659     if (evt_reco_track_gr_All -> GetN() != 0) {evt_reco_track_gr_All -> Set(0);}
1660     if (evt_reco_track_gr_PhiLoose -> GetN() != 0) {evt_reco_track_gr_PhiLoose -> Set(0);}
1661     if (evt_reco_track_gr_Z -> GetN() != 0) {evt_reco_track_gr_Z -> Set(0);}
1662     if (evt_reco_track_gr_ZDCA -> GetN() != 0) {evt_reco_track_gr_ZDCA -> Set(0);}
1663     if (evt_reco_track_gr_ZDCAPhi -> GetN() != 0) {evt_reco_track_gr_ZDCAPhi -> Set(0);}
1664     if (evt_true_track_gr -> GetN() != 0) {evt_true_track_gr -> Set(0);}
1665     if (track_gr -> GetN() != 0) {track_gr -> Set(0);}
1666
1667     out_recotrack_eta_d.clear();
1668     out_recotrack_phi_d.clear();
1669     out_truetrack_eta_d.clear();
1670     out_truetrack_phi_d.clear();
1671     out_track_eta_i.clear();
1672     out_track_delta_phi_d.clear();
1673
1674     effective_total_NClus = 0;
1675
1676     mega_track_finder -> ClearEvt();
1677
1678     // for (int i = 0; i < 360; i++)
1679     // {
1680     //     inner_clu_phi_map[i].clear();
1681     //     outer_clu_phi_map[i].clear();
1682     // }
1683
1684     // std::fill(std::begin(inner_clu_phi_map), std::end(inner_clu_phi_map), std::vector<pair<bool,clu_info>>());
1685     // std::fill(std::begin(outer_clu_phi_map), std::end(outer_clu_phi_map), std::vector<pair<bool,clu_info>>());
1686     inner_clu_phi_map.clear();
1687     outer_clu_phi_map.clear();
1688     inner_clu_phi_map = vector<vector<pair<bool,clu_info>>>(360);
1689     outer_clu_phi_map = vector<vector<pair<bool,clu_info>>>(360);
1690
1691     proto_pair_index.clear();
1692     proto_pair_delta_phi_abs.clear();
1693     proto_pair_delta_phi.clear();
1694     inner_used_clu.clear();
1695     outer_used_clu.clear();
1696
1697     out_N2Clu_track = 0;
1698     out_N3Clu_track = 0;
1699     out_N4Clu_track = 0;
1700
1701     return_tag = 0;
1702     evt_NTrack = 0;
1703     evt_NTrack_MC = 0;
1704     return;
1705 }
1706
1707
1708 void INTTEta::PrintPlots()
1709 {
1710
1711     // note : plots for the mega track finder
1712
1713     NClu3_track_centrality_2D = mega_track_finder -> Get_NClu3_track_centrality_2D();
1714     NClu4_track_centrality_2D = mega_track_finder -> Get_NClu4_track_centrality_2D();
1715     cluster4_track_phi_1D = mega_track_finder -> Get_cluster4_track_phi_1D();
1716     cluster3_track_phi_1D = mega_track_finder -> Get_cluster3_track_phi_1D();
1717     cluster3_inner_track_eta_1D = mega_track_finder -> Get_cluster3_inner_track_eta_1D();
1718     cluster3_inner_track_phi_1D = mega_track_finder -> Get_cluster3_inner_track_phi_1D();
1719     cluster3_outer_track_eta_1D = mega_track_finder -> Get_cluster3_outer_track_eta_1D();
1720     cluster3_outer_track_phi_1D = mega_track_finder -> Get_cluster3_outer_track_phi_1D();
1721     clu4_track_ReducedChi2_1D = mega_track_finder -> Get_clu4_track_ReducedChi2_1D();
1722     clu3_track_ReducedChi2_1D = mega_track_finder -> Get_clu3_track_ReducedChi2_1D();
1723     mega_track_eta_1D = mega_track_finder -> Get_mega_track_eta_1D();
1724
1725     c1 -> cd();
1726
1727     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1728     mega_track_finding_ratio_2D -> Draw("colz0");
1729     mega_track_finding_ratio_2D -> SetMinimum(0);
1730     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1731     c1 -> Print( Form("%s/mega_track_finding_ratio_2D.pdf", out_folder_directory.c_str()) );
1732     c1 -> Clear();
1733
1734
1735     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1736     cluster3_inner_track_eta_1D -> Draw("hist");
1737     cluster3_inner_track_eta_1D -> SetMinimum(0);
1738     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1739     c1 -> Print( Form("%s/cluster3_inner_track_eta_1D.pdf", out_folder_directory.c_str()) );
1740     c1 -> Clear();
1741
1742     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1743     cluster3_outer_track_eta_1D -> Draw("hist");
1744     cluster3_outer_track_eta_1D -> SetMinimum(0);
1745     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1746     c1 -> Print( Form("%s/cluster3_outer_track_eta_1D.pdf", out_folder_directory.c_str()) );
1747     c1 -> Clear();
1748
1749     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1750     mega_track_eta_1D -> Draw("hist");
1751     mega_track_eta_1D -> SetMinimum(0);
1752     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1753     c1 -> Print( Form("%s/mega_track_eta_1D.pdf", out_folder_directory.c_str()) );
1754     c1 -> Clear();
1755
1756     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1757     check_inner_layer_clu_phi_1D -> Draw("hist");
1758     check_inner_layer_clu_phi_1D -> SetMinimum(0);
1759     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1760     c1 -> Print( Form("%s/check_inner_layer_clu_phi_1D.pdf", out_folder_directory.c_str()) );
1761     c1 -> Clear();
1762
1763     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1764     check_outer_layer_clu_phi_1D -> Draw("hist");
1765     check_outer_layer_clu_phi_1D -> SetMinimum(0);
1766     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1767     c1 -> Print( Form("%s/check_outer_layer_clu_phi_1D.pdf", out_folder_directory.c_str()) );
1768     c1 -> Clear();
1769
1770     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1771     clu4_track_ReducedChi2_1D -> Draw("hist");
1772     clu4_track_ReducedChi2_1D -> SetMinimum(0);
1773     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1774     c1 -> Print( Form("%s/clu4_track_ReducedChi2_1D.pdf", out_folder_directory.c_str()) );
1775     c1 -> Clear();
1776
1777     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1778     clu3_track_ReducedChi2_1D -> Draw("hist");
1779     clu3_track_ReducedChi2_1D -> SetMinimum(0);
1780     clu3_track_ReducedChi2_1D -> SetMaximum(clu3_track_ReducedChi2_1D->GetBinContent(clu3_track_ReducedChi2_1D->GetMaximumBin())*1.3);
1781     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1782     coord_line -> DrawLine(mega_track_finder->Get_performance_cut(),0,mega_track_finder->Get_performance_cut(),clu3_track_ReducedChi2_1D->GetBinContent(clu3_track_ReducedChi2_1D->GetMaximumBin())*1.3);
1783     c1 -> Print( Form("%s/clu3_track_ReducedChi2_1D.pdf", out_folder_directory.c_str()) );
1784     c1 -> Clear();
1785
1786     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1787     c1 -> cd();
1788     c1 -> Print( Form("%s/track_cluster_ratio_1D.pdf(", out_folder_directory.c_str()) );
1789     c1 -> Clear();
1790     for (int i = 0; i < track_cluster_ratio_1D.size(); i++)
1791     {
1792         c1 -> cd();
1793         track_cluster_ratio_1D[i] -> Draw("hist");
1794         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1795         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
1796         c1 -> Print(Form("%s/track_cluster_ratio_1D.pdf", out_folder_directory.c_str()));
1797         c1 -> Clear();
1798     }
1799     c1 -> Print( Form("%s/track_cluster_ratio_1D.pdf)", out_folder_directory.c_str()) );
1800
1801     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1802     c1 -> cd();
1803     c1 -> Print( Form("%s/track_cluster_ratio_1D_MC.pdf(", out_folder_directory.c_str()) );
1804     c1 -> Clear();
1805     for (int i = 0; i < track_cluster_ratio_1D_MC.size(); i++)
1806     {
1807         c1 -> cd();
1808         track_cluster_ratio_1D_MC[i] -> Draw("hist");
1809         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1810         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
1811         c1 -> Print(Form("%s/track_cluster_ratio_1D_MC.pdf", out_folder_directory.c_str()));
1812         c1 -> Clear();
1813     }
1814     c1 -> Print( Form("%s/track_cluster_ratio_1D_MC.pdf)", out_folder_directory.c_str()) );
1815
1816     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1817     // c1 -> Print( Form("%s/dNdeta_1D.pdf(", out_folder_directory.c_str()) );
1818     // for (int i = 0; i < dNdeta_1D.size(); i++)
1819     // {
1820     //     double N_correction_evt = (i == dNdeta_1D_MC.size() - 1) ? N_GoodEvent : N_GoodEvent_vec[i];
1821     //     dNdeta_1D_MC[i] -> Scale(1./double(dNdeta_1D_MC[i] -> GetBinWidth(1) ));
1822     //     dNdeta_1D_MC[i] -> Scale(1./double(N_correction_evt));
1823
1824     //     dNdeta_1D[i] -> Scale(1./double(dNdeta_1D[i] -> GetBinWidth(1) ));
1825     //     dNdeta_1D[i] -> Scale(1./double(N_correction_evt));
1826     //     dNdeta_1D[i] -> GetYaxis() -> SetRangeUser(0,800);
1827     //     dNdeta_1D[i] -> Draw("ep");
1828     //     dNdeta_1D_MC[i] -> Draw("hist same");
1829     //     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1830     //     draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s, inclusive Zvtx",centrality_region[i].c_str()));
1831     //     draw_text -> DrawLatex(0.21, 0.86, Form("Nevt MC = Nevt reco: %.1f", N_correction_evt));
1832     //     draw_text -> DrawLatex(0.21, 0.82, Form("MC entry : %.2f, tight entry : %.2f", dNdeta_1D_MC[i] -> Integral(), dNdeta_1D[i] -> Integral()));
1833     //     c1 -> Print(Form("%s/dNdeta_1D.pdf", out_folder_directory.c_str()));
1834     //     c1 -> Clear();
1835     // }
1836     // c1 -> Print( Form("%s/dNdeta_1D.pdf)", out_folder_directory.c_str()) );
1837
1838     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1839     c1 -> Print( Form("%s/track_eta_phi_2D.pdf(", out_folder_directory.c_str()) );
1840     for (int i = 0; i < track_eta_phi_2D.size(); i++)
1841     {
1842         track_eta_phi_2D[i] -> Draw("colz0");
1843         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1844         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
1845         c1 -> Print(Form("%s/track_eta_phi_2D.pdf", out_folder_directory.c_str()));
1846         c1 -> Clear();
1847     }
1848     c1 -> Print( Form("%s/track_eta_phi_2D.pdf)", out_folder_directory.c_str()) );
1849
1850     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1851     c1 -> Print( Form("%s/track_eta_z_2D.pdf(", out_folder_directory.c_str()) );
1852     for (int i = 0; i < track_eta_z_2D.size(); i++)
1853     {
1854         track_eta_z_2D[i] -> Draw("colz0");
1855         DrawEtaZGrid();
1856         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1857         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
1858         c1 -> Print(Form("%s/track_eta_z_2D.pdf", out_folder_directory.c_str()));
1859         c1 -> Clear();
1860     }
1861     c1 -> Print( Form("%s/track_eta_z_2D.pdf)", out_folder_directory.c_str()) );
1862
1863     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1864     c1 -> Print( Form("%s/track_delta_phi_1D.pdf(", out_folder_directory.c_str()) );
1865     for (int i = 0; i < track_delta_phi_1D.size(); i++)
1866     {
1867         track_delta_phi_1D[i] -> SetMinimum(0);
1868         track_delta_phi_1D[i] -> Draw("hist");
1869         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1870         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
1871         coord_line -> DrawLine(-1*phi_diff_cut, 0, -1 * phi_diff_cut, track_delta_phi_1D[i]->GetMaximum() * 1.1);
1872         coord_line -> DrawLine(phi_diff_cut, 0, phi_diff_cut, track_delta_phi_1D[i]->GetMaximum() * 1.1);
1873         c1 -> Print(Form("%s/track_delta_phi_1D.pdf", out_folder_directory.c_str()));
1874         c1 -> Clear();
1875     }
1876     c1 -> Print( Form("%s/track_delta_phi_1D.pdf)", out_folder_directory.c_str()) );
1877
1878     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1879     c1 -> Print( Form("%s/track_delta_eta_1D.pdf(", out_folder_directory.c_str()) );
1880     for (int i = 0; i < track_delta_eta_1D.size(); i++)
1881     {
1882         track_delta_eta_1D[i] -> SetMinimum(0);
1883         track_delta_eta_1D[i] -> Draw("hist");
1884         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1885         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
1886         // coord_line -> DrawLine(-1*phi_diff_cut, 0, -1 * phi_diff_cut, track_delta_eta_1D[i]->GetMaximum() * 1.1);
1887         // coord_line -> DrawLine(phi_diff_cut, 0, phi_diff_cut, track_delta_eta_1D[i]->GetMaximum() * 1.1);
1888         c1 -> Print(Form("%s/track_delta_eta_1D.pdf", out_folder_directory.c_str()));
1889         c1 -> Clear();
1890     }
1891     c1 -> Print( Form("%s/track_delta_eta_1D.pdf)", out_folder_directory.c_str()) );
1892
1893     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1894     c1 -> Print( Form("%s/track_delta_eta_1D_post.pdf(", out_folder_directory.c_str()) );
1895     for (int i = 0; i < track_delta_eta_1D_post.size(); i++)
1896     {
1897         track_delta_eta_1D_post[i] -> SetMinimum(0);
1898         track_delta_eta_1D_post[i] -> Draw("hist");
1899         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1900         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
1901         // coord_line -> DrawLine(-1*phi_diff_cut, 0, -1 * phi_diff_cut, track_delta_eta_1D_post[i]->GetMaximum() * 1.1);
1902         // coord_line -> DrawLine(phi_diff_cut, 0, phi_diff_cut, track_delta_eta_1D_post[i]->GetMaximum() * 1.1);
1903         c1 -> Print(Form("%s/track_delta_eta_1D_post.pdf", out_folder_directory.c_str()));
1904         c1 -> Clear();
1905     }
1906     c1 -> Print( Form("%s/track_delta_eta_1D_post.pdf)", out_folder_directory.c_str()) );
1907
1908     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1909     c1 -> Print( Form("%s/track_DCA_distance.pdf(", out_folder_directory.c_str()) );
1910     for (int i = 0; i < track_DCA_distance.size(); i++)
1911     {
1912         track_DCA_distance[i] -> SetMinimum(0);
1913         track_DCA_distance[i] -> Draw("hist");
1914         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1915         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
1916         coord_line -> DrawLine(DCA_cut.first, 0, DCA_cut.first, track_DCA_distance[i]->GetMaximum() * 1.05);
1917         coord_line -> DrawLine(DCA_cut.second, 0, DCA_cut.second, track_DCA_distance[i]->GetMaximum() * 1.05);
1918         c1 -> Print(Form("%s/track_DCA_distance.pdf", out_folder_directory.c_str()));
1919         c1 -> Clear();
1920     }
1921     c1 -> Print( Form("%s/track_DCA_distance.pdf)", out_folder_directory.c_str()) );
1922
1923     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1924     c1 -> Print( Form("%s/track_phi_DCA_2D.pdf(", out_folder_directory.c_str()) );
1925     for (int i = 0; i < track_phi_DCA_2D.size(); i++)
1926     {
1927         track_phi_DCA_2D[i] -> Draw("colz0");
1928         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1929         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
1930
1931         coord_line -> DrawLine(track_phi_DCA_2D[i] -> GetXaxis() -> GetXmin(), DCA_cut.first, track_phi_DCA_2D[i] -> GetXaxis() -> GetXmax(), DCA_cut.first);
1932         coord_line -> DrawLine(track_phi_DCA_2D[i] -> GetXaxis() -> GetXmin(), DCA_cut.second, track_phi_DCA_2D[i] -> GetXaxis() -> GetXmax(), DCA_cut.second);
1933         coord_line -> DrawLine(-1 * phi_diff_cut, track_phi_DCA_2D[i] -> GetYaxis() -> GetXmin(), -1 * phi_diff_cut, track_phi_DCA_2D[i] -> GetYaxis() -> GetXmax());
1934         coord_line -> DrawLine(phi_diff_cut, track_phi_DCA_2D[i] -> GetYaxis() -> GetXmin(), phi_diff_cut, track_phi_DCA_2D[i] -> GetYaxis() -> GetXmax());
1935         c1 -> Print(Form("%s/track_phi_DCA_2D.pdf", out_folder_directory.c_str()));
1936         c1 -> Clear();
1937     }
1938     c1 -> Print( Form("%s/track_phi_DCA_2D.pdf)", out_folder_directory.c_str()) );
1939
1940     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1941     c1 -> Print( Form("%s/track_ratio_1D.pdf(", out_folder_directory.c_str()) );
1942     for (int i = 0; i < track_ratio_1D.size(); i++)
1943     {
1944         track_ratio_1D[i] -> SetMinimum(0);
1945         track_ratio_1D[i] -> Draw("hist");
1946         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
1947         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
1948         c1 -> Print(Form("%s/track_ratio_1D.pdf", out_folder_directory.c_str()));
1949         c1 -> Clear();
1950     }
1951     c1 -> Print( Form("%s/track_ratio_1D.pdf)", out_folder_directory.c_str()) );
1952
1953     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
1954     for (int i = 0; i < final_track_delta_phi_1D.size(); i++)
1955     {
1956         c1 -> Print( Form("%s/final_track_delta_phi_1D_MBin%i.pdf(", out_folder_directory.c_str(), i) );
1957         for (int i1 = 0; i1 < final_track_delta_phi_1D[i].size(); i1++)
1958         {
1959             if (final_track_delta_phi_1D[i][i1] -> GetEntries() > 0) {
1960                 double hist_offset = get_dist_offset(final_track_delta_phi_1D[i][i1], 15);
1961                 gaus_pol1_fit->SetParameters( final_track_delta_phi_1D[i][i1] -> GetBinContent(final_track_delta_phi_1D[i][i1] -> GetMaximumBin()) - hist_offset, 0, final_track_delta_phi_1D[i][i1]->GetStdDev()/2., hist_offset, 0);
1962                 final_track_delta_phi_1D[i][i1] -> Fit(gaus_pol1_fit,"NQ");
1963
1964
1965                 // note : par[0] : size
1966                 // note : par[1] : ratio of the two gaussians
1967                 // note : par[2] : mean
1968                 // note : par[3] : width of gaus 1
1969                 // note : par[4] : width of gaus 2
1970                 // note : par[5] : offset
1971                 // note : par[6] : slope
1972                 // note : fit with double gaussian function + pol1
1973                 d_gaus_pol1_fit -> SetParameters(final_track_delta_phi_1D[i][i1] -> GetBinContent(final_track_delta_phi_1D[i][i1] -> GetMaximumBin()) - hist_offset, 0.2, 0, final_track_delta_phi_1D[i][i1]->GetStdDev()/2., final_track_delta_phi_1D[i][i1]->GetStdDev()/2., hist_offset, 0);
1974                 d_gaus_pol1_fit -> SetParLimits(1, 0, 0.5);  // note : the first gaussian is  the main distribution, So it should contain more than 50% of the total distribution
1975                 d_gaus_pol1_fit -> SetParLimits(3, 0, 1000); // note : the width of the gaussian should be positive
1976                 d_gaus_pol1_fit -> SetParLimits(4, 0, 1000); // note : the width of the gaussian should be positive
1977                 final_track_delta_phi_1D[i][i1] -> Fit(d_gaus_pol1_fit,"NQ");
1978                 // note : extract the signal region
1979                 draw_d_gaus -> SetParameters(d_gaus_pol1_fit->GetParameter(0), d_gaus_pol1_fit->GetParameter(1), d_gaus_pol1_fit->GetParameter(2), d_gaus_pol1_fit->GetParameter(3), d_gaus_pol1_fit->GetParameter(4));
1980                 // note : extract the part of pol1 background
1981                 draw_pol1_line -> SetParameters(d_gaus_pol1_fit -> GetParameter(5), d_gaus_pol1_fit -> GetParameter(6));
1982
1983                 // note : fit the background region only by the pol2 function
1984                 // note : p[0] + p[1]*(x-p[3])+p[2] * (x-p[3])^2
1985                 bkg_fit_pol2 -> SetParameters(hist_offset, 0, -0.2, 0, signal_region);
1986                 bkg_fit_pol2 -> FixParameter(4, signal_region);
1987                 bkg_fit_pol2 -> SetParLimits(2, -100, 0);
1988                 final_track_delta_phi_1D[i][i1] -> Fit(bkg_fit_pol2,"NQ");
1989                 // note : extract the background region (which includes the signal region also)
1990                 draw_pol2_line -> SetParameters(bkg_fit_pol2 -> GetParameter(0), bkg_fit_pol2 -> GetParameter(1), bkg_fit_pol2 -> GetParameter(2), bkg_fit_pol2 -> GetParameter(3));
1991             }
1992
1993
1994
1995             final_track_delta_phi_1D[i][i1] -> SetMinimum(0);
1996             final_track_delta_phi_1D[i][i1] -> SetMaximum( final_track_delta_phi_1D[i][i1] -> GetBinContent(final_track_delta_phi_1D[i][i1] -> GetMaximumBin()) * 1.5);
1997             final_track_delta_phi_1D[i][i1] -> Draw("hist");
1998             // d_gaus_pol1_fit -> Draw("lsame");
1999             // draw_d_gaus -> Draw("lsame");
2000             // draw_pol1_line -> Draw("lsame");
2001             if (final_track_delta_phi_1D[i][i1] -> GetEntries() > 0) {draw_pol2_line -> Draw("lsame");}
2002
2003             // gaus_pol1_fit -> Draw("lsame");
2004             // draw_gaus_line -> SetParameters(fabs(gaus_pol1_fit -> GetParameter(0)), gaus_pol1_fit -> GetParameter(1), fabs(gaus_pol1_fit -> GetParameter(2)), 0);
2005             // draw_gaus_line -> Draw("lsame");
2006             // draw_pol1_line -> SetParameters(gaus_pol1_fit -> GetParameter(3), gaus_pol1_fit -> GetParameter(4));
2007             // draw_pol1_line -> Draw("lsame");
2008
2009             // cout<<" "<<endl;
2010             // final_eta_entry[i].push_back((draw_gaus_line -> GetParameter(1) - 3 * draw_gaus_line -> GetParameter(2)), draw_gaus_line -> GetParameter(1) + 3 * draw_gaus_line -> GetParameter(2));
2011             // cout<<i<<" "<<i1<<" gaus fit par  : "<<fabs(gaus_pol1_fit -> GetParameter(0))<<" "<<(gaus_pol1_fit -> GetParameter(1))<<" "<<fabs(gaus_pol1_fit -> GetParameter(2))<<endl;
2012             // double gaus_integral = fabs(draw_gaus_line -> Integral( (draw_gaus_line -> GetParameter(1) - 3 * fabs(draw_gaus_line -> GetParameter(2))), draw_gaus_line -> GetParameter(1) + 3 * fabs(draw_gaus_line -> GetParameter(2)) )) / final_track_delta_phi_1D[i][i1] -> GetBinWidth(1);
2013             // cout<<i<<" "<<i1<<" gaus integral : "<< gaus_integral <<endl;
2014
2015             // double d_gaus_integral = fabs(draw_d_gaus -> Integral( -1. * signal_region, signal_region )) / final_track_delta_phi_1D[i][i1] -> GetBinWidth(1);
2016             // cout<<i<<" "<<i1<<" D-gaus integral : "<< d_gaus_integral <<endl;
2017
2018             double pol2_bkg_integral = fabs(draw_pol2_line -> Integral( -1. * signal_region, signal_region )) / final_track_delta_phi_1D[i][i1] -> GetBinWidth(1);
2019             // cout<<i<<" "<<i1<<" pol2_bkg integral: "<<pol2_bkg_integral<<endl;
2020
2021             int binID_X = eta_z_convert_inverse_map[i1].first;
2022             int binID_Y = eta_z_convert_inverse_map[i1].second;
2023             // final_dNdeta_1D[i]->SetBinContent(i1 + 1, good_tracklet_counting[i][i1] - pol2_bkg_integral );
2024             // final_dNdeta_1D[i]->SetBinContent(i1 + 1, good_tracklet_counting[i][i1] ); // note : no background subtraction
2025             if (binID_Y == tight_zvtx_bin) { final_dNdeta_1D[i]->SetBinContent(binID_X, good_tracklet_counting[i][i1] ); } // note : no background subtraction
2026
2027             coord_line -> DrawLine(-1*signal_region, 0, -1 * signal_region, final_track_delta_phi_1D[i][i1] -> GetBinContent(final_track_delta_phi_1D[i][i1] -> GetMaximumBin()) * 1.5);
2028             coord_line -> DrawLine(signal_region, 0, signal_region, final_track_delta_phi_1D[i][i1] -> GetBinContent(final_track_delta_phi_1D[i][i1] -> GetMaximumBin()) * 1.5);
2029
2030             ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2031             draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s, #eta: %.2f ~ %.2f, Z: %.0f ~ %.0f mm",centrality_region[i].c_str(), eta_region_hist -> GetBinCenter(binID_X) - eta_region_hist -> GetBinWidth(binID_X)/2., eta_region_hist -> GetBinCenter(binID_X) + eta_region_hist -> GetBinWidth(binID_X)/2., coarse_eta_z_map -> GetYaxis() -> GetBinCenter(binID_Y) - coarse_eta_z_map -> GetYaxis() -> GetBinWidth(binID_Y)/2., coarse_eta_z_map -> GetYaxis() -> GetBinCenter(binID_Y) + coarse_eta_z_map -> GetYaxis() -> GetBinWidth(binID_Y)/2.));
2032             draw_text -> DrawLatex(0.21, 0.85, Form("MBin: %i, #eta bin: %i, Z bin: %i", i, binID_X, binID_Y));
2033             // draw_text -> DrawLatex(0.21, 0.85, Form("Guassian integral : %.2f", gaus_integral));
2034             // draw_text -> DrawLatex(0.21, 0.80, Form("D-Guassian integral : %.2f", d_gaus_integral));
2035             // draw_text -> DrawLatex(0.21, 0.80, Form("pol2: %.2f + %.2f(x-%.2f) + %.2f(x-%.2f)^{2}", bkg_fit_pol2 -> GetParameter(0), bkg_fit_pol2 -> GetParameter(1), bkg_fit_pol2 -> GetParameter(3), bkg_fit_pol2 -> GetParameter(2), bkg_fit_pol2 -> GetParameter(3)));
2036             c1 -> Print( Form("%s/final_track_delta_phi_1D_MBin%i.pdf", out_folder_directory.c_str(), i) );
2037             c1 -> Clear();
2038         }
2039         c1 -> Print( Form("%s/final_track_delta_phi_1D_MBin%i.pdf)", out_folder_directory.c_str(), i) );
2040     }
2041
2042     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2043     for (int i = 0; i < final_track_multi_delta_phi_1D.size(); i++)
2044     {
2045         c1 -> Print( Form("%s/final_track_multi_delta_phi_1D_MBin%i.pdf(", out_folder_directory.c_str(), i) );
2046         for (int i1 = 0; i1 < final_track_multi_delta_phi_1D[i].size(); i1++)
2047         {
2048             double hist_offset = get_dist_offset(final_track_multi_delta_phi_1D[i][i1], 15);
2049             gaus_pol1_fit->SetParameters( final_track_multi_delta_phi_1D[i][i1] -> GetBinContent(final_track_multi_delta_phi_1D[i][i1] -> GetMaximumBin()) - hist_offset, 0, final_track_multi_delta_phi_1D[i][i1]->GetStdDev()/2., hist_offset, 0);
2050             if (final_track_multi_delta_phi_1D[i][i1] -> GetEntries() > 0) {final_track_multi_delta_phi_1D[i][i1] -> Fit(gaus_pol1_fit,"NQ");}
2051
2052             // note : par[0] : size
2053             // note : par[1] : ratio of the two gaussians
2054             // note : par[2] : mean
2055             // note : par[3] : width of gaus 1
2056             // note : par[4] : width of gaus 2
2057             // note : par[5] : offset
2058             // note : par[6] : slope
2059             // note : fit with double gaussian function + pol1
2060             d_gaus_pol1_fit -> SetParameters(final_track_multi_delta_phi_1D[i][i1] -> GetBinContent(final_track_multi_delta_phi_1D[i][i1] -> GetMaximumBin()) - hist_offset, 0.2, 0, final_track_multi_delta_phi_1D[i][i1]->GetStdDev()/2., final_track_multi_delta_phi_1D[i][i1]->GetStdDev()/2., hist_offset, 0);
2061             d_gaus_pol1_fit -> SetParLimits(1, 0, 0.5);  // note : the first gaussian is  the main distribution, So it should contain more than 50% of the total distribution
2062             d_gaus_pol1_fit -> SetParLimits(3, 0, 1000); // note : the width of the gaussian should be positive
2063             d_gaus_pol1_fit -> SetParLimits(4, 0, 1000); // note : the width of the gaussian should be positive
2064             if (final_track_multi_delta_phi_1D[i][i1] -> GetEntries() > 0) {final_track_multi_delta_phi_1D[i][i1] -> Fit(d_gaus_pol1_fit,"NQ");}
2065             // note : extract the signal region
2066             draw_d_gaus -> SetParameters(d_gaus_pol1_fit->GetParameter(0), d_gaus_pol1_fit->GetParameter(1), d_gaus_pol1_fit->GetParameter(2), d_gaus_pol1_fit->GetParameter(3), d_gaus_pol1_fit->GetParameter(4));
2067             // note : extract the part of pol1 background
2068             draw_pol1_line -> SetParameters(d_gaus_pol1_fit -> GetParameter(5), d_gaus_pol1_fit -> GetParameter(6));
2069
2070             // note : fit the background region only by the pol2 function
2071             // note : p[0] + p[1]*(x-p[3])+p[2] * (x-p[3])^2
2072             bkg_fit_pol2 -> SetParameters(hist_offset, 0, -0.2, 0, signal_region);
2073             bkg_fit_pol2 -> FixParameter(4, signal_region);
2074             bkg_fit_pol2 -> SetParLimits(2, -100, 0);
2075             if (final_track_multi_delta_phi_1D[i][i1] -> GetEntries() > 0) {final_track_multi_delta_phi_1D[i][i1] -> Fit(bkg_fit_pol2,"NQ");}
2076             // note : extract the background region (which includes the signal region also)
2077             draw_pol2_line -> SetParameters(bkg_fit_pol2 -> GetParameter(0), bkg_fit_pol2 -> GetParameter(1), bkg_fit_pol2 -> GetParameter(2), bkg_fit_pol2 -> GetParameter(3));
2078
2079             final_track_multi_delta_phi_1D[i][i1] -> SetMinimum(0);
2080             final_track_multi_delta_phi_1D[i][i1] -> SetMaximum( final_track_multi_delta_phi_1D[i][i1] -> GetBinContent(final_track_multi_delta_phi_1D[i][i1] -> GetMaximumBin()) * 1.5);
2081             final_track_multi_delta_phi_1D[i][i1] -> Draw("hist");
2082             // d_gaus_pol1_fit -> Draw("lsame");
2083             // draw_d_gaus -> Draw("lsame");
2084             // draw_pol1_line -> Draw("lsame");
2085             draw_pol2_line -> Draw("lsame");
2086
2087             // gaus_pol1_fit -> Draw("lsame");
2088             // draw_gaus_line -> SetParameters(fabs(gaus_pol1_fit -> GetParameter(0)), gaus_pol1_fit -> GetParameter(1), fabs(gaus_pol1_fit -> GetParameter(2)), 0);
2089             // draw_gaus_line -> Draw("lsame");
2090             // draw_pol1_line -> SetParameters(gaus_pol1_fit -> GetParameter(3), gaus_pol1_fit -> GetParameter(4));
2091             // draw_pol1_line -> Draw("lsame");
2092
2093             // cout<<" "<<endl;
2094             // final_eta_entry[i].push_back((draw_gaus_line -> GetParameter(1) - 3 * draw_gaus_line -> GetParameter(2)), draw_gaus_line -> GetParameter(1) + 3 * draw_gaus_line -> GetParameter(2));
2095             // cout<<i<<" "<<i1<<" gaus fit par  : "<<fabs(gaus_pol1_fit -> GetParameter(0))<<" "<<(gaus_pol1_fit -> GetParameter(1))<<" "<<fabs(gaus_pol1_fit -> GetParameter(2))<<endl;
2096             double gaus_integral = fabs(draw_gaus_line -> Integral( (draw_gaus_line -> GetParameter(1) - 3 * fabs(draw_gaus_line -> GetParameter(2))), draw_gaus_line -> GetParameter(1) + 3 * fabs(draw_gaus_line -> GetParameter(2)) )) / final_track_multi_delta_phi_1D[i][i1] -> GetBinWidth(1);
2097             // cout<<i<<" "<<i1<<" gaus integral : "<< gaus_integral <<endl;
2098
2099             double d_gaus_integral = fabs(draw_d_gaus -> Integral( -1. * signal_region, signal_region )) / final_track_multi_delta_phi_1D[i][i1] -> GetBinWidth(1);
2100             // cout<<i<<" "<<i1<<" D-gaus integral : "<< d_gaus_integral <<endl;
2101
2102             double pol2_bkg_integral = fabs(draw_pol2_line -> Integral( -1. * signal_region, signal_region )) / final_track_multi_delta_phi_1D[i][i1] -> GetBinWidth(1);
2103             // cout<<i<<" "<<i1<<" pol2_bkg integral: "<<pol2_bkg_integral<<endl;
2104
2105             // final_dNdeta_multi_1D[i]->SetBinContent(i1 + 1, good_tracklet_multi_counting[i][i1] - pol2_bkg_integral );
2106             int binID_X = eta_z_convert_inverse_map[i1].first;
2107             int binID_Y = eta_z_convert_inverse_map[i1].second;
2108             if (binID_Y == tight_zvtx_bin && final_track_multi_delta_phi_1D[i][i1] -> GetEntries() != 0) { final_dNdeta_multi_1D[i]->SetBinContent(binID_X, good_tracklet_multi_counting[i][i1] - pol2_bkg_integral ); } // note : no background subtraction
2109
2110             coord_line -> DrawLine(-1*signal_region, 0, -1 * signal_region, final_track_multi_delta_phi_1D[i][i1] -> GetBinContent(final_track_multi_delta_phi_1D[i][i1] -> GetMaximumBin()) * 1.5);
2111             coord_line -> DrawLine(signal_region, 0, signal_region, final_track_multi_delta_phi_1D[i][i1] -> GetBinContent(final_track_multi_delta_phi_1D[i][i1] -> GetMaximumBin()) * 1.5);
2112
2113             ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2114             draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s, #eta: %.2f ~ %.2f, Z: %.0f ~ %.0f mm",centrality_region[i].c_str(), eta_region_hist -> GetBinCenter(binID_X) - eta_region_hist -> GetBinWidth(binID_X)/2., eta_region_hist -> GetBinCenter(binID_X) + eta_region_hist -> GetBinWidth(binID_X)/2., coarse_eta_z_map -> GetYaxis() -> GetBinCenter(binID_Y) - coarse_eta_z_map -> GetYaxis() -> GetBinWidth(binID_Y)/2., coarse_eta_z_map -> GetYaxis() -> GetBinCenter(binID_Y) + coarse_eta_z_map -> GetYaxis() -> GetBinWidth(binID_Y)/2.));
2115             draw_text -> DrawLatex(0.21, 0.85, Form("MBin: %i, #eta bin: %i, Z bin: %i", i, binID_X, binID_Y));
2116             // draw_text -> DrawLatex(0.21, 0.85, Form("Guassian integral : %.2f", gaus_integral));
2117             // draw_text -> DrawLatex(0.21, 0.80, Form("D-Guassian integral : %.2f", d_gaus_integral));
2118             draw_text -> DrawLatex(0.21, 0.80, Form("pol2: %.2f + %.2f(x-%.2f) + %.2f(x-%.2f)^{2}", bkg_fit_pol2 -> GetParameter(0), bkg_fit_pol2 -> GetParameter(1), bkg_fit_pol2 -> GetParameter(3), bkg_fit_pol2 -> GetParameter(2), bkg_fit_pol2 -> GetParameter(3)));
2119             c1 -> Print( Form("%s/final_track_multi_delta_phi_1D_MBin%i.pdf", out_folder_directory.c_str(), i) );
2120             c1 -> Clear();
2121         }
2122         c1 -> Print( Form("%s/final_track_multi_delta_phi_1D_MBin%i.pdf)", out_folder_directory.c_str(), i) );
2123     }
2124
2125     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2126     c1 -> Print( Form("%s/final_dNdeta_1D.pdf(", out_folder_directory.c_str()) );
2127     for (int i = 0; i < final_dNdeta_1D.size(); i++)
2128     {
2129         double N_correction_evt_MC = (i == final_dNdeta_1D_MC.size() - 1) ? MBin_Z_evt_hist_2D_MC -> GetBinContent(MBin_Z_evt_hist_2D_MC -> GetNbinsX(),tight_zvtx_bin) : MBin_Z_evt_hist_2D_MC -> GetBinContent(i+1,tight_zvtx_bin);
2130         final_dNdeta_1D_MC[i] -> Scale(1./double(final_dNdeta_1D_MC[i] -> GetBinWidth(1) ));
2131         final_dNdeta_1D_MC[i] -> Scale(1./double(N_correction_evt_MC));
2132
2133         double N_correction_evt = (i == final_dNdeta_1D_MC.size() - 1) ? MBin_Z_evt_hist_2D -> GetBinContent(MBin_Z_evt_hist_2D -> GetNbinsX(),tight_zvtx_bin) : MBin_Z_evt_hist_2D -> GetBinContent(i+1,tight_zvtx_bin);
2134         final_dNdeta_1D[i] -> Scale(1./double(final_dNdeta_1D[i] -> GetBinWidth(1) ));
2135         final_dNdeta_1D[i] -> Scale(1./double(N_correction_evt));
2136         final_dNdeta_1D[i] -> GetYaxis() -> SetRangeUser(0, final_dNdeta_1D[i] -> GetMaximum() * 1.5);
2137
2138         final_dNdeta_multi_1D[i] -> Scale(1./double(final_dNdeta_multi_1D[i] -> GetBinWidth(1) ));
2139         final_dNdeta_multi_1D[i] -> Scale(1./double(N_correction_evt));
2140
2141         // final_dNdeta_1D[i] -> SetMinimum(0);
2142         final_dNdeta_1D[i] -> Draw("ep");
2143         final_dNdeta_1D_MC[i] -> Draw("hist same");
2144         final_dNdeta_multi_1D[i] -> Draw("ep same");
2145         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2146         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s, Z: %.2f ~ %.2f",centrality_region[i].c_str(), coarse_eta_z_map -> GetYaxis() -> GetBinCenter(tight_zvtx_bin) - coarse_eta_z_map -> GetYaxis() -> GetBinWidth(tight_zvtx_bin)/2., coarse_eta_z_map -> GetYaxis() -> GetBinCenter(tight_zvtx_bin) + coarse_eta_z_map -> GetYaxis() -> GetBinWidth(tight_zvtx_bin)/2.));
2147         draw_text -> DrawLatex(0.21, 0.86, Form("Nevt MC : %.1f, Nevt reco : %.1f",N_correction_evt_MC, N_correction_evt));
2148         draw_text -> DrawLatex(0.21, 0.82, Form("MC entry : %.2f, tight entry : %.2f, loose entry : %.2f", final_dNdeta_1D_MC[i] -> Integral(), final_dNdeta_1D[i] -> Integral(), final_dNdeta_multi_1D[i] -> Integral()));
2149         c1 -> Print(Form("%s/final_dNdeta_1D.pdf", out_folder_directory.c_str()));
2150         c1 -> Clear();
2151     }
2152     c1 -> Print( Form("%s/final_dNdeta_1D.pdf)", out_folder_directory.c_str()) );
2153
2154     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2155     track_cluster_ratio_multiplicity_2D -> Draw("colz0");
2156     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2157     c1 -> Print( Form("%s/track_cluster_ratio_multiplicity_2D.pdf", out_folder_directory.c_str()) );
2158     c1 -> Clear();
2159
2160     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2161     track_cluster_ratio_multiplicity_2D_MC -> Draw("colz0");
2162     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2163     c1 -> Print( Form("%s/track_cluster_ratio_multiplicity_2D_MC.pdf", out_folder_directory.c_str()) );
2164     c1 -> Clear();
2165
2166     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2167     track_correlation_2D -> Draw("colz0");
2168     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2169     correlation_Line -> Draw("lsame");
2170     c1 -> Print( Form("%s/track_correlation_2D.pdf", out_folder_directory.c_str()) );
2171     c1 -> Clear();
2172
2173     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2174     track_ratio_2D -> Draw("colz0");
2175     coord_line -> DrawLine( track_ratio_2D->GetXaxis()->GetXmin(), 1, track_ratio_2D->GetXaxis()->GetXmax(), 1 );
2176     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2177     correlation_Line -> Draw("lsame");
2178     c1 -> Print( Form("%s/track_ratio_2D.pdf", out_folder_directory.c_str()) );
2179     c1 -> Clear();
2180
2181     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2182     reco_eta_correlation_2D -> Draw("colz0");
2183     correlation_Line -> Draw("lsame");
2184     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2185     c1 -> Print( Form("%s/reco_eta_correlation_2D.pdf", out_folder_directory.c_str()) );
2186     c1 -> Clear();
2187
2188     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2189     reco_eta_diff_reco3P_2D -> Draw("colz0");
2190     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2191     c1 -> Print( Form("%s/reco_eta_diff_reco3P_2D.pdf", out_folder_directory.c_str()) );
2192     c1 -> Clear();
2193
2194     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2195     reco_eta_diff_1D -> Draw("hist");
2196     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2197     c1 -> Print( Form("%s/reco_eta_diff_1D.pdf", out_folder_directory.c_str()) );
2198     c1 -> Clear();
2199
2200     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2201     c1 -> Print( Form("%s/track_DeltaPhi_eta_2D.pdf(", out_folder_directory.c_str()) );
2202     for (int i = 0; i < track_DeltaPhi_eta_2D.size(); i++)
2203     {
2204         track_DeltaPhi_eta_2D[i] -> Draw("colz0");
2205         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2206         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
2207         c1 -> Print(Form("%s/track_DeltaPhi_eta_2D.pdf", out_folder_directory.c_str()));
2208         c1 -> Clear();
2209     }
2210     c1 -> Print( Form("%s/track_DeltaPhi_eta_2D.pdf)", out_folder_directory.c_str()) );
2211
2212     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2213     c1 -> Print( Form("%s/track_DeltaPhi_DeltaEta_2D.pdf(", out_folder_directory.c_str()) );
2214     for (int i = 0; i < track_DeltaPhi_DeltaEta_2D.size(); i++)
2215     {
2216         track_DeltaPhi_DeltaEta_2D[i] -> Draw("colz0");
2217         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2218         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
2219         c1 -> Print(Form("%s/track_DeltaPhi_DeltaEta_2D.pdf", out_folder_directory.c_str()));
2220         c1 -> Clear();
2221     }
2222     c1 -> Print( Form("%s/track_DeltaPhi_DeltaEta_2D.pdf)", out_folder_directory.c_str()) );
2223
2224     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2225     // c1 -> Print( Form("%s/dNdeta_1D_MC.pdf(", out_folder_directory.c_str()) );
2226     // for (int i = 0; i < dNdeta_1D_MC.size(); i++)
2227     // {
2228     //     // dNdeta_1D_MC[i] -> Scale(1./double(dNdeta_1D_MC[i] -> GetBinWidth(1) ));
2229     //     // double N_correction_evt = (i == dNdeta_1D_MC.size() - 1) ? N_GoodEvent : N_GoodEvent_vec[i];
2230     //     // dNdeta_1D_MC[i] -> Scale(1./double(N_GoodEvent));
2231     //     dNdeta_1D_MC[i] -> GetYaxis() -> SetRangeUser(0,800);
2232     //     dNdeta_1D_MC[i] -> Draw("hist");
2233     //     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2234     //     draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
2235     //     c1 -> Print(Form("%s/dNdeta_1D_MC.pdf", out_folder_directory.c_str()));
2236     //     c1 -> Clear();
2237     // }
2238     // c1 -> Print( Form("%s/dNdeta_1D_MC.pdf)", out_folder_directory.c_str()) );
2239
2240
2241     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2242     c1 -> Print( Form("%s/track_eta_z_2D_MC.pdf(", out_folder_directory.c_str()) );
2243     for (int i = 0; i < track_eta_z_2D_MC.size(); i++)
2244     {
2245         track_eta_z_2D_MC[i] -> Draw("colz0");
2246         DrawEtaZGrid();
2247         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2248         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
2249         c1 -> Print(Form("%s/track_eta_z_2D_MC.pdf", out_folder_directory.c_str()));
2250         c1 -> Clear();
2251     }
2252     c1 -> Print( Form("%s/track_eta_z_2D_MC.pdf)", out_folder_directory.c_str()) );
2253
2254     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2255     clu_used_centrality_2D -> Draw("colz0");
2256     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2257     c1 -> Print( Form("%s/clu_used_centrality_2D.pdf", out_folder_directory.c_str()) );
2258     c1 -> Clear();
2259
2260     // note : for the mega cluster study
2261     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2262     cluster4_track_phi_1D -> Draw("hist");
2263     cluster4_track_phi_1D -> SetMinimum(0);
2264     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2265     c1 -> Print( Form("%s/cluster4_track_phi_1D.pdf", out_folder_directory.c_str()) );
2266     c1 -> Clear();
2267
2268     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2269     cluster3_track_phi_1D -> Draw("hist");
2270     cluster3_track_phi_1D -> SetMinimum(0);
2271     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2272     c1 -> Print( Form("%s/cluster3_track_phi_1D.pdf", out_folder_directory.c_str()) );
2273     c1 -> Clear();
2274
2275     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2276     cluster3_inner_track_phi_1D -> Draw("hist");
2277     cluster3_inner_track_phi_1D -> SetMinimum(0);
2278     for (int i = 0; i < cluster3_inner_track_phi_1D -> GetNbinsX(); i++)
2279     {
2280         if (cluster3_inner_track_phi_1D -> GetBinContent(i+1) > 200)
2281         {
2282             cout<<i+1<<" mega Inner track phi bin center: "<<cluster3_inner_track_phi_1D -> GetBinCenter(i+1)<<" bin content: "<<cluster3_inner_track_phi_1D -> GetBinContent(i+1)<<endl;
2283         }
2284     }
2285
2286     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2287     c1 -> Print( Form("%s/cluster3_inner_track_phi_1D.pdf", out_folder_directory.c_str()) );
2288     c1 -> Clear();
2289
2290     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2291     cluster3_outer_track_phi_1D -> Draw("hist");
2292     cluster3_outer_track_phi_1D -> SetMinimum(0);
2293     for (int i = 0; i < cluster3_outer_track_phi_1D -> GetNbinsX(); i++)
2294     {
2295         if (cluster3_outer_track_phi_1D -> GetBinContent(i+1) > 200)
2296         {
2297             cout<<i+1<<" mega Outer track phi bin center: "<<cluster3_outer_track_phi_1D -> GetBinCenter(i+1)<<" bin content: "<<cluster3_outer_track_phi_1D -> GetBinContent(i+1)<<endl;
2298         }
2299     }
2300     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2301     c1 -> Print( Form("%s/cluster3_outer_track_phi_1D.pdf", out_folder_directory.c_str()) );
2302     c1 -> Clear();
2303
2304     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2305     NClu4_track_centrality_2D -> Draw("colz0");
2306     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2307     c1 -> Print( Form("%s/NClu4_track_centrality_2D.pdf", out_folder_directory.c_str()) );
2308     c1 -> Clear();
2309
2310     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2311     NClu3_track_centrality_2D -> Draw("colz0");
2312     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2313     c1 -> Print( Form("%s/NClu3_track_centrality_2D.pdf", out_folder_directory.c_str()) );
2314     c1 -> Clear();
2315
2316     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2317     coarse_eta_z_map -> Draw("colz0");
2318     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2319     c1 -> Print( Form("%s/coarse_eta_z_map.pdf", out_folder_directory.c_str()) );
2320     c1 -> Clear();
2321
2322     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2323     c1 -> Print( Form("%s/coarse_eta_z_2D_MC.pdf(", out_folder_directory.c_str()) );
2324     for (int i = 0; i < coarse_eta_z_2D_MC.size(); i++)
2325     {
2326         coarse_eta_z_2D_MC[i] -> Draw("colz0");
2327         DrawEtaZGrid();
2328         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2329         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
2330         c1 -> Print(Form("%s/coarse_eta_z_2D_MC.pdf", out_folder_directory.c_str()));
2331         c1 -> Clear();
2332     }
2333     c1 -> Print( Form("%s/coarse_eta_z_2D_MC.pdf)", out_folder_directory.c_str()) );
2334
2335     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2336     MBin_Z_evt_hist_2D -> Draw("colz0");
2337     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2338     c1 -> Print( Form("%s/MBin_Z_evt_hist_2D.pdf", out_folder_directory.c_str()) );
2339     c1 -> Clear();
2340
2341     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2342     MBin_Z_evt_hist_2D_MC -> Draw("colz0");
2343     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2344     c1 -> Print( Form("%s/MBin_Z_evt_hist_2D_MC.pdf", out_folder_directory.c_str()) );
2345     c1 -> Clear();
2346
2347     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2348     c1 -> Print( Form("%s/coarse_Reco_SignalNTracklet_Single_eta_z_2D.pdf(", out_folder_directory.c_str()) );
2349     for (int i = 0; i < coarse_Reco_SignalNTracklet_Single_eta_z_2D.size(); i++)
2350     {
2351         coarse_Reco_SignalNTracklet_Single_eta_z_2D[i] -> Draw("colz0");
2352         DrawEtaZGrid();
2353         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2354         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
2355         c1 -> Print(Form("%s/coarse_Reco_SignalNTracklet_Single_eta_z_2D.pdf", out_folder_directory.c_str()));
2356         c1 -> Clear();
2357     }
2358     c1 -> Print( Form("%s/coarse_Reco_SignalNTracklet_Single_eta_z_2D.pdf)", out_folder_directory.c_str()) );
2359
2360     // note : ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2361     c1 -> Print( Form("%s/coarse_Reco_SignalNTracklet_Multi_eta_z_2D.pdf(", out_folder_directory.c_str()) );
2362     for (int i = 0; i < coarse_Reco_SignalNTracklet_Multi_eta_z_2D.size(); i++)
2363     {
2364         coarse_Reco_SignalNTracklet_Multi_eta_z_2D[i] -> Draw("colz0");
2365         DrawEtaZGrid();
2366         ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, Form("#it{#bf{sPHENIX}} %s", plot_text.c_str()));
2367         draw_text -> DrawLatex(0.21, 0.90, Form("Centrality : %s",centrality_region[i].c_str()));
2368         c1 -> Print(Form("%s/coarse_Reco_SignalNTracklet_Multi_eta_z_2D.pdf", out_folder_directory.c_str()));
2369         c1 -> Clear();
2370     }
2371     c1 -> Print( Form("%s/coarse_Reco_SignalNTracklet_Multi_eta_z_2D.pdf)", out_folder_directory.c_str()) );
2372 }
2373
2374 void INTTEta::EndRun()
2375 {
2376     for (int i = 0; i < N_GoodEvent_vec.size(); i++)
2377     {
2378         cout<<"N good evt break down : "<< i <<" "<<N_GoodEvent_vec[i]<<endl;
2379     }
2380
2381     cout<<"N good evt inclusive : "<<N_GoodEvent<<endl;
2382
2383     out_file -> cd();
2384     tree_out -> SetDirectory(out_file);
2385
2386     // note : write the histogram of DeltaPhi for each centrality bin, each eta bin and each z bin
2387     // note : for the one that allows single cluster getting involved in multiple tracks
2388     for (int i = 0; i < final_track_multi_delta_phi_1D.size(); i++)
2389     {
2390         for (int i1 = 0; i1 < final_track_multi_delta_phi_1D[i].size(); i1++)
2391         {
2392             final_track_multi_delta_phi_1D[i][i1] -> Write(Form("Reco_DeltaPhi1D_Multi_MBin%i_Eta%i_Z%i", i, eta_z_convert_inverse_map[i1].first, eta_z_convert_inverse_map[i1].second));
2393         }
2394     }
2395
2396     // note : write the histogram of DeltaPhi for each centrality bin, each eta bin and each z bin
2397     // note : for the ones that only allow single cluster to be used in one track
2398     for (int i = 0; i < final_track_delta_phi_1D.size(); i++)
2399     {
2400         for (int i1 = 0; i1 < final_track_delta_phi_1D[i].size(); i1++)
2401         {
2402             final_track_delta_phi_1D[i][i1] -> Write(Form("Reco_DeltaPhi1D_Single_MBin%i_Eta%i_Z%i", i, eta_z_convert_inverse_map[i1].first, eta_z_convert_inverse_map[i1].second));
2403         }
2404     }
2405
2406     // note : write the histogram to the root file, the one containing the N reco tracks in signal region
2407     for (int i = 0; i < coarse_Reco_SignalNTracklet_Multi_eta_z_2D.size(); i++)
2408     {
2409         coarse_Reco_SignalNTracklet_Multi_eta_z_2D[i] -> Write(Form("Reco_SignalNTracklet_Multi_MBin%i",i));
2410     }
2411     for (int i = 0; i < coarse_Reco_SignalNTracklet_Single_eta_z_2D.size(); i++)
2412     {
2413         coarse_Reco_SignalNTracklet_Single_eta_z_2D[i] -> Write(Form("Reco_SignalNTracklet_Single_MBin%i",i));
2414     }
2415
2416     // note : number of true track for each centrality bin, each eta bin and each z bin
2417     for (int i = 0; i < coarse_eta_z_2D_MC.size(); i++)
2418     {
2419         coarse_eta_z_2D_MC[i] -> Write(Form("NTrueTrack_MBin%i",i));
2420     }
2421
2422     // note :
2423     for (int i = 0; i < coarse_eta_z_2D_fulleta_MC.size(); i++)
2424     {
2425         coarse_eta_z_2D_fulleta_MC[i] -> Write(Form("NTrueTrack_FullEta_MBin%i",i));
2426     }
2427
2428     // note : the map for the eta-z reference
2429     coarse_eta_z_map -> Write("Eta_Z_reference");
2430
2431     // note : the number of event for each centrality bin and z vertex bin
2432     MBin_Z_evt_hist_2D -> Write("MBin_Z_evt_map");
2433
2434     // note : the number of event for each centrality bin and z vertex bin, the zvtx is from MC
2435     MBin_Z_evt_hist_2D_MC -> Write("MBin_Z_evt_map_MC");
2436
2437     // note : to save the fine binning 2D eta-z track historgram
2438     for (int i = 0; i < track_eta_z_2D_MC.size(); i++)
2439     {
2440         track_eta_z_2D_MC[i] -> Write(Form("Fine_TrueTrackEtaZ_MBin%i",i));
2441     }
2442
2443     for (int i = 0; i < track_eta_z_2D.size(); i++)
2444     {
2445         track_eta_z_2D[i] -> Write(Form("Fine_RecoTrackEtaZ_MBin%i",i));
2446     }
2447
2448     // note : to save the 2D histogram, the clu_used_centrality_2D
2449     clu_used_centrality_2D -> Write("CluUsed_Centrality2D");
2450
2451     // note : for the mega tracklet study
2452     // note : save the 2D histogram for
2453     cluster3_track_phi_1D -> Write("cluster3_track_phi_1D");
2454     cluster3_inner_track_phi_1D -> Write("cluster3_inner_track_phi_1D");
2455     cluster3_outer_track_phi_1D -> Write("cluster3_outer_track_phi_1D");
2456     cluster4_track_phi_1D -> Write("cluster4_track_phi_1D");
2457     NClu3_track_centrality_2D -> Write("NClu3_track_centrality_2D");
2458     NClu4_track_centrality_2D -> Write("NClu4_track_centrality_2D");
2459     clu3_track_ReducedChi2_1D -> Write("clu3_track_ReducedChi2_1D");
2460     clu4_track_ReducedChi2_1D -> Write("clu4_track_ReducedChi2_1D");
2461
2462     for (int i = 0; i < track_cluster_ratio_1D_MC.size(); i++){track_cluster_ratio_1D_MC[i] -> Write(Form("track_cluster_ratio_1D_MC_MBin%i",i));}
2463     for (int i = 0; i < track_cluster_ratio_1D.size(); i++) {track_cluster_ratio_1D[i] -> Write(Form("track_cluster_ratio_1D_MBin%i",i));}
2464
2465     track_cluster_ratio_multiplicity_2D_MC -> Write("track_cluster_ratio_multiplicity_2D_MC");
2466     track_cluster_ratio_multiplicity_2D -> Write("track_cluster_ratio_multiplicity_2D");
2467     track_correlation_2D -> Write("track_correlation_2D");
2468
2469     for (int i = 0; i < track_DeltaPhi_eta_2D.size(); i++){ track_DeltaPhi_eta_2D[i] -> Write(Form("track_DeltaPhi_eta_2D_MBin%i",i)); }
2470     for (int i = 0; i < track_eta_phi_2D.size(); i++){ track_eta_phi_2D[i] -> Write(Form("track_eta_phi_2D_MBin%i",i)); }
2471     for (int i = 0; i < track_ratio_1D.size(); i++){ track_ratio_1D[i] -> Write(Form("track_ratio_1D_MBin%i",i)); }
2472     track_ratio_2D -> Write("track_ratio_2D");
2473
2474     for (int i = 0; i < dNdeta_1D.size(); i++) {dNdeta_1D[i] -> Write(Form("FineBin_Ntracklet_MBin%i",i));}
2475     for (int i = 0; i < dNdeta_1D_MC.size(); i++) {dNdeta_1D_MC[i] -> Write(Form("FineBin_NTrueTrack_MBin%i",i));}
2476     for (int i = 0; i < dNdeta_1D_MC_edge_eta_cut.size(); i++) {dNdeta_1D_MC_edge_eta_cut[i] -> Write(Form("FineBin_NTrueTrack_EdgeEtaCut_MBin%i",i));}
2477
2478     exclusive_NClus_inner -> Write("exclusive_NClus_inner");
2479     exclusive_NClus_outer -> Write("exclusive_NClus_outer");
2480     exclusive_NClus_sum -> Write("exclusive_NClus_sum");
2481
2482     exclusive_cluster_inner_eta -> Write("exclusive_cluster_inner_eta");
2483     exclusive_cluster_inner_phi -> Write("exclusive_cluster_inner_phi");
2484     exclusive_cluster_outer_eta -> Write("exclusive_cluster_outer_eta");
2485     exclusive_cluster_outer_phi -> Write("exclusive_cluster_outer_phi");
2486     exclusive_cluster_all_eta -> Write("exclusive_cluster_all_eta");
2487     exclusive_cluster_all_phi -> Write("exclusive_cluster_all_phi");
2488
2489     exclusive_cluster_inner_eta_phi_2D -> Write("exclusive_cluster_inner_eta_phi_2D");
2490     exclusive_cluster_outer_eta_phi_2D -> Write("exclusive_cluster_outer_eta_phi_2D");
2491
2492     exclusive_tight_tracklet_eta -> Write("exclusive_tight_tracklet_eta");
2493     exclusive_tight_tracklet_phi -> Write("exclusive_tight_tracklet_phi");
2494     exclusive_loose_tracklet_eta -> Write("exclusive_loose_tracklet_eta");
2495     exclusive_loose_tracklet_phi -> Write("exclusive_loose_tracklet_phi");
2496
2497     exclusive_cluster_inner_eta_adc_2D -> Write("exclusive_cluster_inner_eta_adc_2D");
2498     exclusive_cluster_outer_eta_adc_2D -> Write("exclusive_cluster_outer_eta_adc_2D");
2499
2500     exclusive_cluster_inner_adc -> Write("exclusive_cluster_inner_adc");
2501     exclusive_cluster_outer_adc -> Write("exclusive_cluster_outer_adc");
2502
2503     tree_out -> Write("", TObject::kOverwrite);
2504
2505     out_file -> Close();
2506
2507     std::cout<<"in N event : "<<N_GoodEvent<<"the total N reco tracks is short in : "<<N_reco_cluster_short<<std::endl;
2508
2509     return;
2510 }
2511
2512 //note : accumulate the number of entries from both sides of the histogram
2513 // double INTTEta::get_dist_offset(TH1F * hist_in, int check_N_bin) // note : check_N_bin 1 to N bins of hist
2514 // {
2515 //     if (check_N_bin < 0 || check_N_bin > hist_in -> GetNbinsX()) {cout<<" wrong check_N_bin "<<endl; exit(1);}
2516 //     double total_entry = 0;
2517 //     for (int i = 0; i < check_N_bin; i++)
2518 //     {
2519 //         total_entry += hist_in -> GetBinContent(i+1); // note : 1, 2, 3.....
2520 //         total_entry += hist_in -> GetBinContent(hist_in -> GetNbinsX() - i);
2521 //     }
2522
2523 //     return total_entry/double(2. * check_N_bin);
2524 // }
2525
2526
2527 double INTTEta::grEY_stddev(TGraphErrors * input_grr)
2528 {
2529     vector<double> input_vector; input_vector.clear();
2530     for (int i = 0; i < input_grr -> GetN(); i++) {input_vector.push_back( input_grr -> GetPointY(i) );}
2531
2532     double average = accumulate( input_vector.begin(), input_vector.end(), 0.0 ) / double(input_vector.size());
2533
2534     double sum_subt = 0;
2535     for (int ele : input_vector) {sum_subt+=pow((ele-average),2);}
2536
2537     return sqrt(sum_subt/(input_vector.size()-1));
2538 }
2539
2540
2541 pair<double, double> INTTEta::mirrorPolynomial(double a, double b) {
2542     // Interchange 'x' and 'y'
2543     double mirroredA = 1.0 / a;
2544     double mirroredB = -b / a;
2545
2546     return {mirroredA, mirroredB};
2547 }
2548
2549 // note : pair<reduced chi2, eta of the track>
2550 // note : vector : {r,z}
2551 // note : p0 : vertex point {r,z,zerror}
2552 // note : p1 : inner layer
2553 // note : p2 : outer layer
2554 pair<double,double> INTTEta::Get_eta(vector<double>p0, vector<double>p1, vector<double>p2)
2555 {
2556     // if (track_gr -> GetN() != 0) {cout<<"In INTTEta class, track_gr is not empty, track_gr size : "<<track_gr -> GetN()<<endl; exit(0);}
2557
2558     if (track_gr -> GetN() != 0) {track_gr -> Set(0);}
2559
2560     vector<double> r_vec  = {p0[0], p1[0], p2[0]};
2561     vector<double> z_vec  = {p0[1], p1[1], p2[1]};
2562     vector<double> re_vec = {0, 0, 0};
2563     vector<double> ze_vec = {p0[2], ( fabs( p1[1] ) < 130 ) ? 8. : 10., ( fabs( p2[1] ) < 130 ) ? 8. : 10.};
2564
2565     // note : swap the r and z, to have easier fitting
2566     // note : in principle, Z should be in the x axis, R should be in the Y axis
2567     for (int i = 0; i < 3; i++)
2568     {
2569         track_gr -> SetPoint(i,r_vec[i],z_vec[i]);
2570         track_gr -> SetPointError(i,re_vec[i],ze_vec[i]);
2571
2572         // cout<<"In INTTEta class, point : "<<i<<" r : "<<r_vec[i]<<" z : "<<z_vec[i]<<" re : "<<re_vec[i]<<" ze : "<<ze_vec[i]<<endl;
2573     }
2574
2575     double vertical_line = ( fabs( grEY_stddev(track_gr) ) < 0.00001 ) ? 1 : 0;
2576
2577     if (vertical_line) {return {0,0};}
2578     else
2579     {
2580         fit_rz -> SetParameters(0,0);
2581
2582         track_gr -> Fit(fit_rz,"NQ");
2583
2584         pair<double,double> ax_b = mirrorPolynomial(fit_rz -> GetParameter(1),fit_rz -> GetParameter(0));
2585
2586         return  {(fit_rz -> GetChisquare() / double(fit_rz -> GetNDF())), -1 * TMath::Log( fabs( tan( atan2(ax_b.first, (ax_b.first > 0) ? 1. : -1. ) / 2 ) ) )};
2587     }
2588
2589 }
2590
2591 // note : angle_1 = inner clu phi
2592 // note: angle_2 = outer clu phi
2593 double INTTEta::get_delta_phi(double angle_1, double angle_2)
2594 {
2595     vector<double> vec_abs = {fabs(angle_1 - angle_2), fabs(angle_1 - angle_2 + 360), fabs(angle_1 - angle_2 - 360)};
2596     vector<double> vec = {(angle_1 - angle_2), (angle_1 - angle_2 + 360), (angle_1 - angle_2 - 360)};
2597     return vec[std::distance(vec_abs.begin(), std::min_element(vec_abs.begin(),vec_abs.end()))];
2598 }
2599
2600 double INTTEta::get_track_phi(double inner_clu_phi_in, double delta_phi_in)
2601 {
2602     double track_phi = inner_clu_phi_in - (delta_phi_in/2.);
2603     if (track_phi < 0) {track_phi += 360;}
2604     else if (track_phi > 360) {track_phi -= 360;}
2605     else if (track_phi == 360) {track_phi = 0;}
2606     else {track_phi = track_phi;}
2607     return track_phi;
2608 }
2609
2610 double INTTEta::convertTo360(double radian) {
2611
2612     double angle = radian * (180./M_PI);
2613
2614     if (fabs(radian) == M_PI) {return 180.;}
2615     else if (radian == 0)
2616     {
2617         return 0;
2618     }
2619     else if ( radian > 0 )
2620     {
2621         return angle;
2622     }
2623     else
2624     {
2625         return angle + 360;
2626     }
2627 }
2628
2629 double INTTEta::get_clu_eta(vector<double> vertex, vector<double> clu_pos)
2630 {
2631     double correct_x = clu_pos[0] - vertex[0];
2632     double correct_y = clu_pos[1] - vertex[1];
2633     double correct_z = clu_pos[2] - vertex[2];
2634     double clu_r = sqrt(pow(correct_x,2) + pow(correct_y,2));
2635     // cout<<"correct info : "<<correct_x<<" "<<correct_y<<" "<<correct_z<<" "<<clu_r<<endl;
2636
2637     return -0.5 * TMath::Log((sqrt(pow(correct_z,2)+pow(clu_r,2))-(correct_z)) / (sqrt(pow(correct_z,2)+pow(clu_r,2))+(correct_z)));
2638 }
2639
2640 pair<int,int> INTTEta::GetTH2BinXY(int histNbinsX, int histNbinsY, int binglobal)
2641 {
2642     // cout<<"so detail 1 "<<endl;
2643     if (binglobal == -1) {return {-1,-1};}
2644     // cout<<"so detail 2 "<<endl;
2645     int nx  = histNbinsX+2;
2646     int ny  = histNbinsY+2;
2647     // cout<<"so detail 3 "<<endl;
2648     int binx = binglobal%nx;
2649     int biny = ((binglobal-binx)/nx)%ny;
2650     // cout<<"so detail 4 "<<endl;
2651     // cout<<binx<<" "<<biny<<endl;
2652
2653     // note : the binx and biny start from 1
2654     return {binx, biny};
2655 }
2656
2657 double INTTEta::GetTH2Index1D(pair<int,int> XY_index, int histNbinsX)
2658 {
2659     if (XY_index.first == -1) {return -1;}
2660
2661     // note : XY_index.first : x index, starts from 1
2662     // note : XY_index.second : y index, starts from 1
2663     return (XY_index.first - 1) + (XY_index.second - 1) * histNbinsX;
2664 }
2665
2666 void INTTEta::DrawEtaZGrid()
2667 {
2668     // note : draw the vertical line, to segment the eta region
2669     for (int i = 0; i < eta_region.size(); i++) { coord_line -> DrawLine(eta_region[i], z_region[0], eta_region[i], z_region[z_region.size() - 1]); }
2670
2671     // note : draw the horizontal line, to segment the z region
2672     for (int i = 0; i < z_region.size(); i++) { coord_line -> DrawLine(eta_region[0], z_region[i], eta_region[eta_region.size() - 1], z_region[i]); }
2673 }
2674
2675 double INTTEta::get_TH1F_Entries(TH1F * hist_in)
2676 {
2677     double entry_counting = 0;
2678     for (int i = 0; i < hist_in -> GetNbinsX(); i++) {entry_counting += hist_in -> GetBinContent(i+1);}
2679
2680     return entry_counting;
2681 }
2682
2683 #endif
2684
2685
2686 // // note : this function is for the event by event vertex calculation
2687 // // note : this function is the new method, which means that we first put the cluster into a phi map, and then there are two for loops still, but we only check the certain range of the phi
2688 // // note : try the method that allows single cluster to be used multiple times
2689 // // note : so the background subtraction is needed
2690 // void INTTEta::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, int centrality_bin, vector<vector<float>> true_track_info){ // note : evt_z : {z, width}
2691 //     return_tag = 0;
2692
2693 //     if (event_i%1 == 0) {cout<<"In INTTEta class, running event : "<<event_i<<endl;}
2694
2695 //     inner_NClu = temp_sPH_inner_nocolumn_vec.size();
2696 //     outer_NClu = temp_sPH_outer_nocolumn_vec.size();
2697 //     total_NClus = inner_NClu + outer_NClu;
2698
2699 //     // cout<<"test_0"<<endl;
2700 //     if (total_NClus < zvtx_cal_require) {return; cout<<"return confirmation"<<endl;}
2701
2702 //     if (run_type == "MC" && NvtxMC != 1) { return; cout<<"In INTTEta class, event : "<<event_i<<" Nvtx : "<<NvtxMC<<" Nvtx more than one "<<endl;}
2703 //     if (PhiCheckTag == false)            { return; cout<<"In INTTEta class, event : "<<event_i<<" Nvtx : "<<NvtxMC<<" Not full phi has hits "<<endl;}
2704
2705 //     if (inner_NClu < 10 || outer_NClu < 10 || total_NClus > N_clu_cut || total_NClus < N_clu_cutl)
2706 //     {
2707 //         return;
2708 //         printf("In INTTEta class, event : %i, low clu continue, NClus : %lu \n", event_i, total_NClus);
2709 //     }
2710
2711 //     // todo : the z vertex range is here
2712 //     if (-220 > evt_z.first + evt_z.second || -180 < evt_z.first - evt_z.second) {return;}
2713 //     // if (-100 > evt_z.first + evt_z.second || 100 < evt_z.first - evt_z.second) {return;}
2714
2715
2716 //     N_GoodEvent += 1;
2717 //     N_GoodEvent_vec[centrality_map[centrality_bin]] += 1;
2718
2719 //     // cout<<"N inner cluster : "<<inner_NClu<<" N outer cluster : "<<outer_NClu<<endl;
2720
2721 //     double INTT_eta_acceptance_l = -0.5 * TMath::Log((sqrt(pow(-230.-TrigvtxMC[2]*10.,2)+pow(INTT_layer_R[3],2))-(-230.-TrigvtxMC[2]*10.)) / (sqrt(pow(-230.-TrigvtxMC[2]*10.,2)+pow(INTT_layer_R[3],2))+(-230.-TrigvtxMC[2]*10.))); // note : left
2722 //     double INTT_eta_acceptance_r =  -0.5 * TMath::Log((sqrt(pow(230.-TrigvtxMC[2]*10.,2)+pow(INTT_layer_R[3],2))-(230.-TrigvtxMC[2]*10.)) / (sqrt(pow(230.-TrigvtxMC[2]*10.,2)+pow(INTT_layer_R[3],2))+(230.-TrigvtxMC[2]*10.))); // note : right
2723
2724 //     if (run_type == "MC")
2725 //     {
2726 //         // note : for the true track case ----------------------------------------------------------------------------------------------------------------------------------------------------------------
2727 //         // if (event_i % 100 == 0){cout<<"z : "<<TrigvtxMC[2]*10.<<" eta : "<<INTT_eta_acceptance_l<<" "<<INTT_eta_acceptance_r<<endl;}
2728
2729 //         // cout<<"true_track_info : "<<true_track_info.size()<<endl;
2730 //         for (int track_i = 0; track_i < true_track_info.size(); track_i++)
2731 //         {
2732 //             if (true_track_info[track_i][2] == 111 || true_track_info[track_i][2] == 22 || abs(true_track_info[track_i][2]) == 2112){continue;}
2733
2734 //             if (true_track_info[track_i][0] > INTT_eta_acceptance_l && true_track_info[track_i][0] < INTT_eta_acceptance_r)
2735 //             {
2736 //                 dNdeta_1D_MC[centrality_map[centrality_bin]] -> Fill(true_track_info[track_i][0]);
2737 //                 dNdeta_1D_MC[dNdeta_1D_MC.size() - 1]        -> Fill(true_track_info[track_i][0]);
2738 //                 final_dNdeta_1D_MC[centrality_map[centrality_bin]] -> Fill(true_track_info[track_i][0]);
2739 //                 final_dNdeta_1D_MC[dNdeta_1D_MC.size() - 1]        -> Fill(true_track_info[track_i][0]);
2740 //                 track_eta_z_2D_MC[centrality_map[centrality_bin]] -> Fill(true_track_info[track_i][0], TrigvtxMC[2]*10.);
2741 //                 track_eta_z_2D_MC[track_eta_z_2D_MC.size() - 1]   -> Fill(true_track_info[track_i][0], TrigvtxMC[2]*10.);
2742
2743 //                 // cout<<"true track eta : "<<true_track_info[track_i][0]<<" phi : "<<convertTo360(true_track_info[track_i][1])<<endl;
2744 //                 // cout<<"("<<true_track_info[track_i][0]<<", "<<convertTo360(true_track_info[track_i][1])<<"), ";
2745
2746 //                 evt_true_track_gr -> SetPoint(evt_true_track_gr->GetN(),true_track_info[track_i][0], convertTo360(true_track_info[track_i][1]));
2747
2748 //                 evt_NTrack_MC += 1;
2749 //             }
2750 //         }
2751 //         // if (evt_NTrack_MC < 10) {cout<<"evt : "<<event_i<<" ---- N reco track : "<<evt_NTrack<<" N true track : "<<evt_NTrack_MC<<" ratio : "<<double(evt_NTrack) / double(evt_NTrack_MC)<<endl;}
2752 //     }
2753
2754 //     // note : put the cluster into the phi map, the first bool is for the cluster usage.
2755 //     // note : false means the cluster is not used
2756 //     for (int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++) {
2757 //         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());
2758 //         inner_clu_phi_map[ int(Clus_InnerPhi_Offset) ].push_back({false,temp_sPH_inner_nocolumn_vec[inner_i]});
2759
2760 //         double clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{temp_sPH_inner_nocolumn_vec[inner_i].x, temp_sPH_inner_nocolumn_vec[inner_i].y, temp_sPH_inner_nocolumn_vec[inner_i].z});
2761 //         if (clu_eta > INTT_eta_acceptance_l && clu_eta < INTT_eta_acceptance_r) {effective_total_NClus += 1;}
2762 //     }
2763 //     for (int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++) {
2764 //         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());
2765 //         outer_clu_phi_map[ int(Clus_OuterPhi_Offset) ].push_back({false,temp_sPH_outer_nocolumn_vec[outer_i]});
2766
2767 //         double clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{temp_sPH_outer_nocolumn_vec[outer_i].x, temp_sPH_outer_nocolumn_vec[outer_i].y, temp_sPH_outer_nocolumn_vec[outer_i].z});
2768 //         if (clu_eta > INTT_eta_acceptance_l && clu_eta < INTT_eta_acceptance_r) {effective_total_NClus += 1;}
2769 //     }
2770
2771 //     // // note : find the Mega cluster preparation for inner barrel
2772 //     // for (int inner_phi_i = 0; inner_phi_i < 360; inner_phi_i++)
2773 //     // {
2774 //     //     for (int inner_phi_clu_i = 0; inner_phi_clu_i < inner_clu_phi_map[inner_phi_i].size(); inner_phi_clu_i++)
2775 //     //     {
2776 //     //         // note : if the cluster is used, then skip
2777 //     //         if (inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].first == true) {continue;}
2778
2779 //     //         Clus_InnerPhi_Offset_1 = (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());
2780
2781 //     //         // note : the scan range of the -1, 0, 1
2782 //     //         for (int scan_i = -1; scan_i < 2; scan_i++)
2783 //     //         {
2784 //     //             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;
2785 //     //             for (int inner_scan_clu_i = 0; inner_scan_clu_i < inner_clu_phi_map[true_scan_i].size(); inner_scan_clu_i++)
2786 //     //             {
2787 //     //                 // note : if the cluster is used, then skip
2788 //     //                 if (inner_clu_phi_map[true_scan_i][inner_scan_clu_i].first == true) {continue;}
2789
2790 //     //                 // note : the cluster itself
2791 //     //                 if (true_scan_i == inner_phi_i && inner_phi_clu_i == inner_scan_clu_i) {continue;}
2792
2793 //     //                 // note : if it has the same sub layer ID, then skip
2794 //     //                 // todo : this cut may only work for the MC, for the data, it requires additional check
2795 //     //                 if (inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.layer == inner_clu_phi_map[true_scan_i][inner_scan_clu_i].second.layer) {continue;}
2796
2797 //     //                 // note : I expect the two cluster have to have the same z position
2798 //     //                 // note : but in case of data, the z position of the same strip may be fluctuated a little bit
2799 //     //                 // todo : currently, set the Z position flutuation to be 4 mm
2800 //     //                 if (fabs(inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z - inner_clu_phi_map[true_scan_i][inner_scan_clu_i].second.z) > 4 ) {continue;}
2801
2802 //     //                 Clus_InnerPhi_Offset_2 = (inner_clu_phi_map[true_scan_i][inner_scan_clu_i].second.y - beam_origin.second < 0) ? atan2(inner_clu_phi_map[true_scan_i][inner_scan_clu_i].second.y - beam_origin.second, inner_clu_phi_map[true_scan_i][inner_scan_clu_i].second.x - beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(inner_clu_phi_map[true_scan_i][inner_scan_clu_i].second.y - beam_origin.second, inner_clu_phi_map[true_scan_i][inner_scan_clu_i].second.x - beam_origin.first) * (180./TMath::Pi());
2803
2804 //     //                 double delta_phi = get_delta_phi(Clus_InnerPhi_Offset_1, Clus_InnerPhi_Offset_2);
2805
2806 //     //                 // note : if the two cluster are too far away in phi angle, the skip
2807 //     //                 // todo : the delta phi cut may need to be adjusted
2808 //     //                 if (fabs( delta_phi ) > 0.6) {continue;}
2809
2810 //     //                 // note : the pair that passed the cut mentioned above can be considered as a proto_mega_cluster
2811 //     //                 Mega_inner_clu_pair_index.push_back({{inner_phi_i, inner_phi_clu_i, true_scan_i, inner_scan_clu_i}});
2812 //     //                 Mega_inner_clu_delta_phi_abs.push_back(fabs(delta_phi));
2813 //     //                 Mega_inner_clu_phi.push_back( get_track_phi(Clus_InnerPhi_Offset_1, delta_phi) );
2814 //     //             }
2815 //     //         }
2816 //     //     }
2817 //     // }
2818
2819 //     // // note : finr the Mega cluste preparation for outer barrel
2820 //     // for (int outer_phi_i = 0; outer_phi_i < 360; outer_phi_i++)
2821 //     // {
2822 //     //     for (int outer_phi_clu_i = 0; outer_phi_clu_i < outer_clu_phi_map[outer_phi_i].size(); outer_phi_clu_i++)
2823 //     //     {
2824 //     //         // note : if the cluster is used, then skip
2825 //     //         if (outer_clu_phi_map[outer_phi_i][outer_phi_clu_i].first == true) {continue;}
2826
2827 //     //         Clus_OuterPhi_Offset_1 = (outer_clu_phi_map[outer_phi_i][outer_phi_clu_i].second.y - beam_origin.second < 0) ? atan2(outer_clu_phi_map[outer_phi_i][outer_phi_clu_i].second.y - beam_origin.second, outer_clu_phi_map[outer_phi_i][outer_phi_clu_i].second.x - beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(outer_clu_phi_map[outer_phi_i][outer_phi_clu_i].second.y - beam_origin.second, outer_clu_phi_map[outer_phi_i][outer_phi_clu_i].second.x - beam_origin.first) * (180./TMath::Pi());
2828
2829 //     //         // note : the scan range of the -1, 0, 1
2830 //     //         for (int scan_i = -1; scan_i < 2; scan_i++)
2831 //     //         {
2832 //     //             int true_scan_i = ((outer_phi_i + scan_i) < 0) ? 360 + (outer_phi_i + scan_i) : ((outer_phi_i + scan_i) > 359) ? (outer_phi_i + scan_i)-360 : outer_phi_i + scan_i;
2833 //     //             for (int outer_scan_clu_i = 0; outer_scan_clu_i < outer_clu_phi_map[true_scan_i].size(); outer_scan_clu_i++)
2834 //     //             {
2835 //     //                 // note : if the cluster is used, then skip
2836 //     //                 if (outer_clu_phi_map[true_scan_i][outer_scan_clu_i].first == true) {continue;}
2837
2838 //     //                 // note : the cluster itself
2839 //     //                 if (true_scan_i == outer_phi_i && outer_phi_clu_i == outer_scan_clu_i) {continue;}
2840
2841 //     //                 // note : if it has the same sub layer ID, then skip
2842 //     //                 // todo : this cut may only work for the MC, for the data, it requires additional check
2843 //     //                 if (outer_clu_phi_map[outer_phi_i][outer_phi_clu_i].second.layer == outer_clu_phi_map[true_scan_i][outer_scan_clu_i].second.layer) {continue;}
2844
2845 //     //                 // note : I expect the two cluster have to have the same z position
2846 //     //                 // note : but in case of data, the z position of the same strip may be fluctuated a little bit
2847 //     //                 // todo : currently, set the Z position flutuation to be 4 mm
2848 //     //                 if (fabs(outer_clu_phi_map[outer_phi_i][outer_phi_clu_i].second.z - outer_clu_phi_map[true_scan_i][outer_scan_clu_i].second.z) > 4 ) {continue;}
2849
2850 //     //                 Clus_OuterPhi_Offset_2 = (outer_clu_phi_map[true_scan_i][outer_scan_clu_i].second.y - beam_origin.second < 0) ? atan2(outer_clu_phi_map[true_scan_i][outer_scan_clu_i].second.y - beam_origin.second, outer_clu_phi_map[true_scan_i][outer_scan_clu_i].second.x - beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(outer_clu_phi_map[true_scan_i][outer_scan_clu_i].second.y - beam_origin.second, outer_clu_phi_map[true_scan_i][outer_scan_clu_i].second.x - beam_origin.first) * (180./TMath::Pi());
2851
2852 //     //                 double delta_phi = get_delta_phi(Clus_OuterPhi_Offset_1, Clus_OuterPhi_Offset_2);
2853
2854 //     //                 // note : if the two cluster are too far away in phi angle, the skip
2855 //     //                 // todo : the delta phi cut may need to be adjusted
2856 //     //                 if (fabs(delta_phi) > 0.6) {continue;}
2857
2858 //     //                 // note : the pair that passed the cut mentioned above can be considered as a proto_mega_cluster
2859 //     //                 Mega_outer_clu_pair_index.push_back({{outer_phi_i, outer_phi_clu_i, true_scan_i, outer_scan_clu_i}});
2860 //     //                 Mega_outer_clu_delta_phi_abs.push_back(fabs(delta_phi));
2861 //     //                 Mega_outer_clu_phi.push_back( get_track_phi(Clus_OuterPhi_Offset_1, delta_phi) );
2862 //     //             }
2863 //     //         }
2864 //     //     }
2865 //     // }
2866
2867 //     // // note : try to sort the pair of the Inner Mega cluster by the abs delta phi
2868 //     // long long inner_mega_vec_size = Mega_inner_clu_delta_phi_abs.size();
2869 //     // long long ind_inner_mega[Mega_inner_clu_delta_phi_abs.size()];
2870 //     // TMath::Sort(inner_mega_vec_size, &Mega_inner_clu_delta_phi_abs[0], ind_inner_mega, false);
2871
2872 //     // // note : try to sort the pair of the Outer Mega cluster by the abs delta phi
2873 //     // long long outer_mega_vec_size = Mega_outer_clu_delta_phi_abs.size();
2874 //     // long long ind_outer_mega[Mega_outer_clu_delta_phi_abs.size()];
2875 //     // TMath::Sort(outer_mega_vec_size, &Mega_outer_clu_delta_phi_abs[0], ind_outer_mega, false);
2876
2877 //     // // note : 4 cluster track, Inner mega - Outer mega
2878 //     // for (int inner_i; inner_i < Mega_inner_clu_pair_index.size(); inner_i++)
2879 //     // {
2880 //     //     int id_inner[4] = {
2881 //     //         Mega_inner_clu_pair_index[ind_inner_mega[inner_i]][0],
2882 //     //         Mega_inner_clu_pair_index[ind_inner_mega[inner_i]][1], // note : one inner cluster
2883 //     //         Mega_inner_clu_pair_index[ind_inner_mega[inner_i]][2],
2884 //     //         Mega_inner_clu_pair_index[ind_inner_mega[inner_i]][3] // note : another inner cluster
2885 //     //     };
2886
2887 //     //     if (inner_used_mega_clu[Form("%i_%i", id_inner[0], id_inner[1])] != 0) {continue;}
2888 //     //     if (inner_used_mega_clu[Form("%i_%i", id_inner[2], id_inner[3])] != 0) {continue;}
2889
2890 //     //     for (int outer_i; outer_i < Mega_outer_clu_pair_index.size(); outer_i++)
2891 //     //     {
2892 //     //         int id_outer[4] = {
2893 //     //             Mega_outer_clu_pair_index[ind_outer_mega[outer_i]][0],
2894 //     //             Mega_outer_clu_pair_index[ind_outer_mega[outer_i]][1], //  note : one outer cluster
2895 //     //             Mega_outer_clu_pair_index[ind_outer_mega[outer_i]][2],
2896 //     //             Mega_outer_clu_pair_index[ind_outer_mega[outer_i]][3] // note : abother outer cluster
2897 //     //         };
2898
2899 //     //         if (outer_used_mega_clu[Form("%i_%i", id_outer[0], id_outer[1])] != 0) {continue;}
2900 //     //         if (outer_used_mega_clu[Form("%i_%i", id_outer[2], id_outer[3])] != 0) {continue;}
2901
2902 //     //         // note : the delta phi cut
2903 //     //         // todo : the delta phi cut value is fixed here
2904 //     //         if (fabs(get_delta_phi(Mega_inner_clu_phi[ind_inner_mega[inner_i]], Mega_outer_clu_phi[ind_outer_mega[outer_i]])) > 0.6) {continue;}
2905
2906 //     //         int id_inner_small_r = ( get_radius(inner_clu_phi_map[id_inner[0]][id_inner[1]].second.x - beam_origin.first, inner_clu_phi_map[id_inner[0]][id_inner[1]].second.y - beam_origin.second) < get_radius(inner_clu_phi_map[id_inner[2]][id_inner[3]].second.x - beam_origin.first, inner_clu_phi_map[id_inner[2]][id_inner[3]].second.y - beam_origin.second) ) ? 0 : 2;
2907 //     //         int id_outer_big_r   = ( get_radius(outer_clu_phi_map[id_outer[0]][id_outer[1]].second.x - beam_origin.first, outer_clu_phi_map[id_outer[0]][id_outer[1]].second.y - beam_origin.second) > get_radius(outer_clu_phi_map[id_outer[2]][id_outer[3]].second.x - beam_origin.first, outer_clu_phi_map[id_outer[2]][id_outer[3]].second.y - beam_origin.second) ) ? 0 : 2;
2908
2909 //     //         pair<double,double> z_range_info = Get_possible_zvtx(
2910 //     //             0., // get_radius(beam_origin.first,beam_origin.second),
2911 //     //             {get_radius(inner_clu_phi_map[id_inner[id_inner_small_r]][id_inner[id_inner_small_r+1]].second.x - beam_origin.first, inner_clu_phi_map[id_inner[id_inner_small_r]][id_inner[id_inner_small_r+1]].second.y - beam_origin.second), inner_clu_phi_map[id_inner[id_inner_small_r]][id_inner[id_inner_small_r+1]].second.z}, // note : unsign radius
2912 //     //             {get_radius(outer_clu_phi_map[id_outer[id_outer_big_r]][id_outer[id_outer_big_r+1]].second.x - beam_origin.first, outer_clu_phi_map[id_outer[id_outer_big_r]][id_outer[id_outer_big_r+1]].second.y - beam_origin.second), outer_clu_phi_map[id_outer[id_outer_big_r]][id_outer[id_outer_big_r+1]].second.z}  // note : unsign radius
2913 //     //         );
2914
2915 //     //         // 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
2916 //     //         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;}
2917
2918 //     //         cluster4_track_phi_1D -> Fill(get_track_phi(Mega_inner_clu_phi[ind_inner_mega[inner_i]], get_delta_phi(Mega_inner_clu_phi[ind_inner_mega[inner_i]], Mega_outer_clu_phi[ind_outer_mega[outer_i]])));
2919
2920 //     //         NClu4_track_count += 1;
2921
2922 //     //         // note : to mark the cluster as used
2923 //     //         inner_used_mega_clu[Form("%i_%i", id_inner[0], id_inner[1])] += 1;
2924 //     //         inner_used_mega_clu[Form("%i_%i", id_inner[2], id_inner[3])] += 1;
2925 //     //         outer_used_mega_clu[Form("%i_%i", id_outer[0], id_outer[1])] += 1;
2926 //     //         outer_used_mega_clu[Form("%i_%i", id_outer[2], id_outer[3])] += 1;
2927
2928 //     //         inner_clu_phi_map[id_inner[0]][id_inner[1]].first = true;
2929 //     //         inner_clu_phi_map[id_inner[2]][id_inner[3]].first = true;
2930 //     //         outer_clu_phi_map[id_outer[0]][id_outer[1]].first = true;
2931 //     //         outer_clu_phi_map[id_outer[2]][id_outer[3]].first = true;
2932 //     //     }
2933 //     // }
2934
2935 //     // // note : 3 cluster track, Inner mega - Outer normal
2936 //     // for (int inner_i; inner_i < Mega_inner_clu_pair_index.size(); inner_i++) // note : the inner mega cluster loop
2937 //     // {
2938 //     //     int id_inner[4] = {
2939 //     //         Mega_inner_clu_pair_index[ind_inner_mega[inner_i]][0],
2940 //     //         Mega_inner_clu_pair_index[ind_inner_mega[inner_i]][1], // note : one inner cluster
2941 //     //         Mega_inner_clu_pair_index[ind_inner_mega[inner_i]][2],
2942 //     //         Mega_inner_clu_pair_index[ind_inner_mega[inner_i]][3] // note : another inner cluster
2943 //     //     };
2944
2945 //     //     if (inner_used_mega_clu[Form("%i_%i", id_inner[0], id_inner[1])] != 0) {continue;}
2946 //     //     if (inner_used_mega_clu[Form("%i_%i", id_inner[2], id_inner[3])] != 0) {continue;}
2947
2948 //     //     if (Mega_inner_clu_phi[ind_inner_mega[inner_i]] < 0 || Mega_inner_clu_phi[ind_inner_mega[inner_i]] >= 360) {cout<<"test: the mega inner cluster phi angle calculation is wrong, the value is : "<<Mega_inner_clu_phi[ind_inner_mega[inner_i]]<<endl;} // todo : this is a debug line
2949 //     //     int mega_clu_phi_index = int(Mega_inner_clu_phi[ind_inner_mega[inner_i]]);
2950
2951 //     //     // note : outer cluster loop -1, 0, 1
2952 //     //     for (int scan_i = -1; scan_i < 2; scan_i++)
2953 //     //     {
2954 //     //         int true_scan_i = ((mega_clu_phi_index + scan_i) < 0) ? 360 + (mega_clu_phi_index + scan_i) : ((mega_clu_phi_index + scan_i) > 359) ? (mega_clu_phi_index + scan_i)-360 : mega_clu_phi_index + scan_i;
2955
2956 //     //         for (int outer_i; outer_i < outer_clu_phi_map[true_scan_i].size(); outer_i++)
2957 //     //         {
2958 //     //             // note : if the cluster is used, then skip
2959 //     //             if (outer_clu_phi_map[true_scan_i][outer_i].first == true) {continue;}
2960
2961 //     //             // note : calculate the outer cluster phi with the consideration of the beam position
2962 //     //             Clus_OuterPhi_Offset = (outer_clu_phi_map[true_scan_i][outer_i].second.y - beam_origin.second < 0) ? atan2(outer_clu_phi_map[true_scan_i][outer_i].second.y - beam_origin.second, outer_clu_phi_map[true_scan_i][outer_i].second.x - beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(outer_clu_phi_map[true_scan_i][outer_i].second.y - beam_origin.second, outer_clu_phi_map[true_scan_i][outer_i].second.x - beam_origin.first) * (180./TMath::Pi());
2963
2964 //     //             // note : the delta phi cut between the inner-mega-cluster and the outer-cluster
2965 //     //             if ( fabs(get_delta_phi(Mega_inner_clu_phi[ind_inner_mega[inner_i]], Clus_OuterPhi_Offset)) > 0.6 ) { continue; }
2966
2967 //     //             int id_inner_small_r = ( get_radius(inner_clu_phi_map[id_inner[0]][id_inner[1]].second.x - beam_origin.first, inner_clu_phi_map[id_inner[0]][id_inner[1]].second.y - beam_origin.second) < get_radius(inner_clu_phi_map[id_inner[2]][id_inner[3]].second.x - beam_origin.first, inner_clu_phi_map[id_inner[2]][id_inner[3]].second.y - beam_origin.second) ) ? 0 : 2;
2968
2969 //     //             pair<double,double> z_range_info = Get_possible_zvtx(
2970 //     //                 0., // get_radius(beam_origin.first,beam_origin.second),
2971 //     //                 {get_radius(inner_clu_phi_map[id_inner[id_inner_small_r]][id_inner[id_inner_small_r+1]].second.x - beam_origin.first, inner_clu_phi_map[id_inner[id_inner_small_r]][id_inner[id_inner_small_r+1]].second.y - beam_origin.second), inner_clu_phi_map[id_inner[id_inner_small_r]][id_inner[id_inner_small_r+1]].second.z}, // note : unsign radius
2972 //     //                 {get_radius(outer_clu_phi_map[true_scan_i][outer_i].second.x - beam_origin.first, outer_clu_phi_map[true_scan_i][outer_i].second.y - beam_origin.second), outer_clu_phi_map[true_scan_i][outer_i].second.z}  // note : unsign radius
2973 //     //             );
2974
2975 //     //             // 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
2976 //     //             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;}
2977
2978 //     //             cluster3_all_track_phi_1D   -> Fill(get_track_phi(Mega_inner_clu_phi[ind_inner_mega[inner_i]], get_delta_phi(Mega_inner_clu_phi[ind_inner_mega[inner_i]], Clus_OuterPhi_Offset)));
2979 //     //             cluster3_inner_track_phi_1D -> Fill(get_track_phi(Mega_inner_clu_phi[ind_inner_mega[inner_i]], get_delta_phi(Mega_inner_clu_phi[ind_inner_mega[inner_i]], Clus_OuterPhi_Offset)));
2980
2981 //     //             NClu3_track_count += 1;
2982
2983 //     //             // note : to mark the cluster as used
2984 //     //             inner_used_mega_clu[Form("%i_%i", id_inner[0], id_inner[1])] += 1;
2985 //     //             inner_used_mega_clu[Form("%i_%i", id_inner[2], id_inner[3])] += 1;
2986
2987 //     //             inner_clu_phi_map[id_inner[0]][id_inner[1]].first = true;
2988 //     //             inner_clu_phi_map[id_inner[2]][id_inner[3]].first = true;
2989 //     //             outer_clu_phi_map[true_scan_i][outer_i].first = true;
2990
2991 //     //         }
2992 //     //     }
2993 //     // }
2994
2995 //     // // note : 3 cluster track, Inner normal - Outer mega
2996 //     // for (int outer_i; outer_i < Mega_outer_clu_pair_index.size(); outer_i++) // note : the outer mega cluster loop
2997 //     // {
2998 //     //     int id_outer[4] = {
2999 //     //         Mega_outer_clu_pair_index[ind_outer_mega[outer_i]][0],
3000 //     //         Mega_outer_clu_pair_index[ind_outer_mega[outer_i]][1], // note : one outer cluster
3001 //     //         Mega_outer_clu_pair_index[ind_outer_mega[outer_i]][2],
3002 //     //         Mega_outer_clu_pair_index[ind_outer_mega[outer_i]][3] // note : another outer cluster
3003 //     //     };
3004
3005 //     //     if (outer_used_mega_clu[Form("%i_%i", id_outer[0], id_outer[1])] != 0) {continue;}
3006 //     //     if (outer_used_mega_clu[Form("%i_%i", id_outer[2], id_outer[3])] != 0) {continue;}
3007
3008 //     //     if (Mega_outer_clu_phi[ind_outer_mega[outer_i]] < 0 || Mega_outer_clu_phi[ind_outer_mega[outer_i]] >= 360) {cout<<"test: the mega outer cluster phi angle calculation is wrong, the value is : "<<Mega_outer_clu_phi[ind_outer_mega[outer_i]]<<endl;} // todo : this is a debug line
3009 //     //     int mega_clu_phi_index = int(Mega_outer_clu_phi[ind_outer_mega[outer_i]]);
3010
3011 //     //     for (int scan_i = -1; scan_i < 2; scan_i++)
3012 //     //     {
3013 //     //         int true_scan_i = ((mega_clu_phi_index + scan_i) < 0) ? 360 + (mega_clu_phi_index + scan_i) : ((mega_clu_phi_index + scan_i) > 359) ? (mega_clu_phi_index + scan_i)-360 : mega_clu_phi_index + scan_i;
3014
3015 //     //         for (int inner_i; inner_i < inner_clu_phi_map[true_scan_i].size(); inner_i++)
3016 //     //         {
3017 //     //             // note : if the cluster is used, then skip
3018 //     //             if (inner_clu_phi_map[true_scan_i][inner_i].first == true) {continue;}
3019
3020 //     //             // note : calculate the inner cluster phi with the consideration of the beam position
3021 //     //             Clus_InnerPhi_Offset = (inner_clu_phi_map[true_scan_i][inner_i].second.y - beam_origin.second < 0) ? atan2(inner_clu_phi_map[true_scan_i][inner_i].second.y - beam_origin.second, inner_clu_phi_map[true_scan_i][inner_i].second.x - beam_origin.first) * (180./TMath::Pi()) + 360 : atan2(inner_clu_phi_map[true_scan_i][inner_i].second.y - beam_origin.second, inner_clu_phi_map[true_scan_i][inner_i].second.x - beam_origin.first) * (180./TMath::Pi());
3022
3023 //     //             // note : the delta phi cut between the outer-mega-cluster and the inner-cluster
3024 //     //             if ( fabs(get_delta_phi(Mega_outer_clu_phi[ind_outer_mega[outer_i]], Clus_InnerPhi_Offset)) > 0.6 ) { continue; }
3025
3026 //     //             int id_outer_big_r   = ( get_radius(outer_clu_phi_map[id_outer[0]][id_outer[1]].second.x - beam_origin.first, outer_clu_phi_map[id_outer[0]][id_outer[1]].second.y - beam_origin.second) > get_radius(outer_clu_phi_map[id_outer[2]][id_outer[3]].second.x - beam_origin.first, outer_clu_phi_map[id_outer[2]][id_outer[3]].second.y - beam_origin.second) ) ? 0 : 2;
3027
3028 //     //             pair<double,double> z_range_info = Get_possible_zvtx(
3029 //     //                 0., // get_radius(beam_origin.first,beam_origin.second),
3030 //     //                 {get_radius(inner_clu_phi_map[true_scan_i][inner_i].second.x - beam_origin.first, inner_clu_phi_map[true_scan_i][inner_i].second.y - beam_origin.second), inner_clu_phi_map[true_scan_i][inner_i].second.z}, // note : unsign radius
3031 //     //                 {get_radius(outer_clu_phi_map[id_outer[id_outer_big_r]][id_outer[id_outer_big_r+1]].second.x - beam_origin.first, outer_clu_phi_map[id_outer[id_outer_big_r]][id_outer[id_outer_big_r+1]].second.y - beam_origin.second), outer_clu_phi_map[id_outer[id_outer_big_r]][id_outer[id_outer_big_r+1]].second.z}  // note : unsign radius
3032 //     //             );
3033
3034 //     //             // 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
3035 //     //             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;}
3036
3037 //     //             cluster3_all_track_phi_1D   -> Fill(get_track_phi(Clus_InnerPhi_Offset, get_delta_phi(Clus_InnerPhi_Offset, Mega_outer_clu_phi[ind_outer_mega[outer_i]])));
3038 //     //             cluster3_outer_track_phi_1D -> Fill(get_track_phi(Clus_InnerPhi_Offset, get_delta_phi(Clus_InnerPhi_Offset, Mega_outer_clu_phi[ind_outer_mega[outer_i]])));
3039
3040 //     //             NClu3_track_count += 1;
3041
3042 //     //             // note : to mark the cluster as used
3043 //     //             outer_used_mega_clu[Form("%i_%i", id_outer[0], id_outer[1])] += 1;
3044 //     //             outer_used_mega_clu[Form("%i_%i", id_outer[2], id_outer[3])] += 1;
3045
3046 //     //             outer_clu_phi_map[id_outer[0]][id_outer[1]].first = true;
3047 //     //             outer_clu_phi_map[id_outer[2]][id_outer[3]].first = true;
3048 //     //             inner_clu_phi_map[true_scan_i][inner_i].first = true;
3049 //     //         }
3050 //     //     }
3051 //     // }
3052
3053 //     // NClu4_track_centrality_2D -> Fill(centrality_map[centrality_bin], NClu4_track_count);
3054 //     // NClu3_track_centrality_2D -> Fill(centrality_map[centrality_bin], NClu3_track_count);
3055
3056 //     // // note : for the mega cluster test only
3057 //     // return;
3058
3059 //     // note : for two-cluster tracklets only
3060 //     for (int inner_phi_i = 0; inner_phi_i < 360; inner_phi_i++) // note : each phi cell (1 degree)
3061 //     {
3062 //         // note : N cluster in this phi cell
3063 //         for (int inner_phi_clu_i = 0; inner_phi_clu_i < inner_clu_phi_map[inner_phi_i].size(); inner_phi_clu_i++)
3064 //         {
3065 //             clu_multi_used_tight = 0;
3066 //             clu_multi_used_loose = 0;
3067 //             if (inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].first == true) {continue;}
3068
3069 //             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());
3070
3071 //             // todo: change the outer phi scan range
3072 //             // note : the outer phi index, -4, -3, -2, -1, 0, 1, 2, 3, 4
3073 //             for (int scan_i = -4; scan_i < 5; scan_i++)
3074 //             {
3075 //                 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;
3076
3077 //                 // note : N clusters in that outer phi cell
3078 //                 for (int outer_phi_clu_i = 0; outer_phi_clu_i < outer_clu_phi_map[true_scan_i].size(); outer_phi_clu_i++)
3079 //                 {
3080 //                     if (outer_clu_phi_map[true_scan_i][outer_phi_clu_i].first == true) {continue;}
3081
3082 //                     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());
3083 //                     double delta_phi = get_delta_phi(Clus_InnerPhi_Offset, Clus_OuterPhi_Offset);
3084
3085 //                     // if (fabs(delta_phi) > 5.72) {continue;}
3086 //                     if (fabs(delta_phi) > 3.5) {continue;}
3087
3088 //                     double inner_clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{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, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z});
3089 //                     double outer_clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{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, outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.z});
3090 //                     double delta_eta = inner_clu_eta - outer_clu_eta;
3091
3092 //                     track_delta_eta_1D[centrality_map[centrality_bin]] -> Fill( delta_eta );
3093 //                     track_delta_eta_1D[track_delta_eta_1D.size() - 1]  -> Fill( delta_eta );
3094
3095 //                     track_DeltaPhi_DeltaEta_2D[centrality_map[centrality_bin]]        -> Fill(delta_phi, delta_eta);
3096 //                     track_DeltaPhi_DeltaEta_2D[track_DeltaPhi_DeltaEta_2D.size() - 1] -> Fill(delta_phi, delta_eta);
3097
3098 //                     pair<double,double> z_range_info = Get_possible_zvtx(
3099 //                         0., // get_radius(beam_origin.first,beam_origin.second),
3100 //                         {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
3101 //                         {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
3102 //                     );
3103
3104 //                     // 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
3105 //                     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;}
3106 //                     if (fabs(delta_phi) < 0.6) {clu_multi_used_tight += 1;}
3107 //                     track_delta_eta_1D_post[centrality_map[centrality_bin]]     -> Fill( delta_eta );
3108 //                     track_delta_eta_1D_post[track_delta_eta_1D_post.size() - 1] -> Fill( delta_eta );
3109
3110 //                     double DCA_sign = calculateAngleBetweenVectors(
3111 //                         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,
3112 //                         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,
3113 //                         beam_origin.first, beam_origin.second
3114 //                     );
3115
3116 //                     track_DCA_distance[centrality_map[centrality_bin]] -> Fill( DCA_sign );
3117 //                     track_DCA_distance[track_DCA_distance.size() - 1]  -> Fill( DCA_sign );
3118
3119 //                     track_phi_DCA_2D[centrality_map[centrality_bin]] -> Fill( delta_phi, DCA_sign );
3120 //                     track_phi_DCA_2D[track_phi_DCA_2D.size() - 1]    -> Fill( delta_phi, DCA_sign );
3121
3122 //                     track_delta_phi_1D[centrality_map[centrality_bin]] -> Fill( delta_phi );
3123 //                     track_delta_phi_1D[track_delta_phi_1D.size() - 1]  -> Fill( delta_phi );
3124
3125 //                     proto_pair_index.push_back({inner_phi_i, inner_phi_clu_i, true_scan_i, outer_phi_clu_i});
3126 //                     proto_pair_delta_phi_abs.push_back(fabs(delta_phi));
3127 //                     proto_pair_delta_phi.push_back({Clus_InnerPhi_Offset, Clus_OuterPhi_Offset, delta_phi});
3128
3129
3130 //                     // note : do the fill here (find the best match outer cluster with the inner cluster )
3131 //                     Get_eta_pair = Get_eta(
3132 //                         {0., evt_z.first,evt_z.second},
3133 //                         {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},
3134 //                         {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}
3135 //                     );
3136
3137 //                     if  (Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2. > 0.3)
3138 //                     {
3139 //                         cout<<" "<<endl;
3140 //                         cout<<"inner clu eta : "<<inner_clu_eta<<" outer clu eta : "<<outer_clu_eta<<"avg eta : "<< (inner_clu_eta + outer_clu_eta)/2. <<" reco eta : "<<Get_eta_pair.second<<" diff: "<<Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.<<endl;
3141 //                         cout<<"inner clu pos : "<<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<<" "<<inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z<<endl;
3142 //                         cout<<"outer clu pos : "<<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<<" "<<outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second.z<<endl;
3143 //                     }
3144
3145
3146 //                     double track_phi = get_track_phi(Clus_InnerPhi_Offset, delta_phi);
3147
3148 //                     reco_eta_correlation_2D -> Fill(Get_eta_pair.second, (inner_clu_eta + outer_clu_eta)/2.);
3149 //                     reco_eta_diff_reco3P_2D -> Fill(Get_eta_pair.second, Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.);
3150 //                     reco_eta_diff_1D        -> Fill(Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.);
3151
3152 //                     // track_delta_phi_1D[centrality_map[centrality_bin]] -> Fill( delta_phi );
3153 //                     // track_delta_phi_1D[track_delta_phi_1D.size() - 1]  -> Fill( delta_phi );
3154
3155 //                     track_DeltaPhi_eta_2D[centrality_map[centrality_bin]]   -> Fill(delta_phi, Get_eta_pair.second);
3156 //                     track_DeltaPhi_eta_2D[track_DeltaPhi_eta_2D.size() - 1] -> Fill(delta_phi, Get_eta_pair.second);
3157
3158 //                     if (fabs(delta_phi) <= signal_region)
3159 //                     {
3160 //                         dNdeta_1D[centrality_map[centrality_bin]] -> Fill(Get_eta_pair.second);
3161 //                         dNdeta_1D[dNdeta_1D.size() - 1]           -> Fill(Get_eta_pair.second);
3162
3163 //                         track_eta_z_2D[centrality_map[centrality_bin]] -> Fill(Get_eta_pair.second, evt_z.first);
3164 //                         track_eta_z_2D[track_eta_z_2D.size() - 1]      -> Fill(Get_eta_pair.second, evt_z.first);
3165
3166 //                         evt_NTrack += 1;
3167 //                     }
3168
3169 //                     track_eta_phi_2D[centrality_map[centrality_bin]] -> Fill(track_phi, Get_eta_pair.second);
3170 //                     track_eta_phi_2D[track_eta_phi_2D.size() - 1]    -> Fill(track_phi, Get_eta_pair.second);
3171
3172 //                     // note : find the bin of the eta
3173 //                     double eta_bin = eta_region_hist -> Fill(Get_eta_pair.second);
3174 //                     // cout<<"test1 "<<eta_bin<<endl;
3175 //                     if (eta_bin != -1)
3176 //                     {
3177 //                         final_track_delta_phi_1D[centrality_map[centrality_bin]][eta_bin - 1]      -> Fill(delta_phi);
3178 //                         final_track_delta_phi_1D[final_track_delta_phi_1D.size() - 1][eta_bin - 1] -> Fill(delta_phi);
3179
3180 //                         // todo : the signal region selection
3181 //                         if (fabs(delta_phi) <= signal_region) {
3182 //                             good_tracklet_counting[centrality_map[centrality_bin]][eta_bin - 1] += 1;
3183 //                             good_tracklet_counting[final_track_delta_phi_1D.size() - 1][eta_bin - 1] += 1;
3184 //                         }
3185
3186 //                         out_track_delta_phi_d.push_back(delta_phi);
3187 //                         out_recotrack_eta_d.push_back(Get_eta_pair.second);
3188 //                         out_track_eta_i.push_back(eta_bin);
3189 //                     }
3190 //                 }
3191 //             } // note : end outer loop
3192
3193 //             clu_used_centrality_2D -> Fill(total_NClus, clu_multi_used_tight);
3194
3195 //             // // note : if there is no good match with inner and outer clusters, continue
3196 //             // if (clu_multi_used_loose == 0) {continue;}
3197
3198 //             // // note : do the fill here (find the best match outer cluster with the inner cluster )
3199 //             // Get_eta_pair = Get_eta(
3200 //             //     {0., evt_z.first,evt_z.second},
3201 //             //     {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},
3202 //             //     {get_radius(outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.x - beam_origin.first, outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.y - beam_origin.second), outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.z}
3203 //             // );
3204
3205 //             // double inner_clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{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, inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z});
3206 //             // double outer_clu_eta = get_clu_eta({beam_origin.first, beam_origin.second, evt_z.first},{outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.x, outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.y, outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.z});
3207 //             // double track_phi = Clus_InnerPhi_Offset - (pair_delta_phi/2.);
3208
3209 //             // if  (Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2. > 0.3)
3210 //             // {
3211 //             //     cout<<" "<<endl;
3212 //             //     cout<<"inner clu eta : "<<inner_clu_eta<<" outer clu eta : "<<outer_clu_eta<<"avg eta : "<< (inner_clu_eta + outer_clu_eta)/2. <<" reco eta : "<<Get_eta_pair.second<<" diff: "<<Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.<<endl;
3213 //             //     cout<<"inner clu pos : "<<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<<" "<<inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second.z<<endl;
3214 //             //     cout<<"outer clu pos : "<<outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.x<<" "<<outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.y<<" "<<outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].second.z<<endl;
3215 //             // }
3216
3217 //             // reco_eta_correlation_2D -> Fill(Get_eta_pair.second, (inner_clu_eta + outer_clu_eta)/2.);
3218 //             // reco_eta_diff_reco3P_2D -> Fill(Get_eta_pair.second, Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.);
3219 //             // reco_eta_diff_1D        -> Fill(Get_eta_pair.second - (inner_clu_eta + outer_clu_eta)/2.);
3220
3221 //             // // track_delta_phi_1D[centrality_map[centrality_bin]] -> Fill( pair_delta_phi );
3222 //             // // track_delta_phi_1D[track_delta_phi_1D.size() - 1]  -> Fill( pair_delta_phi );
3223
3224 //             // track_DeltaPhi_eta_2D[centrality_map[centrality_bin]]   -> Fill(pair_delta_phi, Get_eta_pair.second);
3225 //             // track_DeltaPhi_eta_2D[track_DeltaPhi_eta_2D.size() - 1] -> Fill(pair_delta_phi, Get_eta_pair.second);
3226
3227 //             // // cout<<"test_5"<<endl;
3228 //             // dNdeta_1D[centrality_map[centrality_bin]] -> Fill(Get_eta_pair.second);
3229 //             // dNdeta_1D[dNdeta_1D.size() - 1]           -> Fill(Get_eta_pair.second);
3230
3231 //             // track_eta_phi_2D[centrality_map[centrality_bin]] -> Fill(track_phi, Get_eta_pair.second);
3232 //             // track_eta_phi_2D[track_eta_phi_2D.size() - 1]    -> Fill(track_phi, Get_eta_pair.second);
3233
3234 //             // track_eta_z_2D[centrality_map[centrality_bin]] -> Fill(Get_eta_pair.second, evt_z.first);
3235 //             // track_eta_z_2D[track_eta_z_2D.size() - 1]      -> Fill(Get_eta_pair.second, evt_z.first);
3236
3237 //             // // note : find the bin of the eta
3238 //             // double eta_bin = eta_region_hist -> Fill(Get_eta_pair.second);
3239 //             // // cout<<"test1 "<<eta_bin<<endl;
3240 //             // if (eta_bin != -1)
3241 //             // {
3242 //             //     final_track_delta_phi_1D[centrality_map[centrality_bin]][eta_bin - 1]      -> Fill(pair_delta_phi);
3243 //             //     final_track_delta_phi_1D[final_track_delta_phi_1D.size() - 1][eta_bin - 1] -> Fill(pair_delta_phi);
3244
3245 //             //     out_track_delta_phi_d.push_back(pair_delta_phi);
3246 //             //     out_recotrack_eta_d.push_back(Get_eta_pair.second);
3247 //             //     out_track_eta_i.push_back(eta_bin);
3248 //             // }
3249
3250 //             // evt_NTrack += 1;
3251 //             // // note : this outer cluster is used
3252 //             // outer_clu_phi_map[pair_outer_index.first][pair_outer_index.second].first = true;
3253 //         }
3254 //     } // note : end inner loop
3255
3256 //     out_eID = event_i;
3257 //     out_evt_centrality_bin = centrality_bin;
3258 //     out_evt_zvtx = evt_z.first;
3259 //     tree_out -> Fill();
3260
3261 //     // cout<<" "<<endl;
3262 //     // cout<<" "<<endl;
3263 //     // cout<<"test_8"<<endl;
3264 //     if (run_type == "MC")
3265 //     {
3266 //         track_correlation_2D -> Fill(evt_NTrack_MC, evt_NTrack);
3267 //         track_ratio_1D[centrality_map[centrality_bin]] -> Fill( double(evt_NTrack) / double(evt_NTrack_MC) );
3268 //         track_ratio_1D[track_ratio_1D.size() - 1]      -> Fill( double(evt_NTrack) / double(evt_NTrack_MC) );
3269 //     }
3270
3271 //     track_cluster_ratio_multiplicity_2D -> Fill( effective_total_NClus, double(effective_total_NClus) / double(evt_NTrack) );
3272 //     track_cluster_ratio_1D[centrality_map[centrality_bin]]    -> Fill( double(effective_total_NClus) / double(evt_NTrack) );
3273 //     track_cluster_ratio_1D[track_cluster_ratio_1D.size() - 1] -> Fill( double(effective_total_NClus) / double(evt_NTrack) );
3274
3275 //     // if (run_type == "MC" && evt_NTrack_MC < draw_evt_cut && event_i % print_plot_ratio == 0){print_evt_plot(event_i, evt_NTrack_MC, inner_NClu, outer_NClu);}
3276 //     return_tag = 1;
3277 // }