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