Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:17:41 

0001 #include "INTTXYvtx.h"
0002 
0003 #include <TF1.h>
0004 
0005 #include <boost/format.hpp>
0006 
0007 #include <cmath>
0008 
0009 double cos_func(double* x, double* par)
0010 {
0011   return -1 * par[0] * cos(par[1] * (x[0] - par[2])) + par[3];
0012 }
0013 
0014 INTTXYvtx::INTTXYvtx(const std::string& runType,
0015                      const std::string& outFolderDirectory,
0016                      std::pair<double, double> beamOrigin,
0017                      double phiDiffCut,
0018                      std::pair<double, double> DCACut,
0019                      int NCluCutl,
0020                      int NCluCut,
0021                      double angleDiffNew_l,
0022                      double angleDiffNew_r,
0023                      double peekCut,
0024                      bool printMessageOpt)
0025   : run_type(runType)
0026   , out_folder_directory(outFolderDirectory)
0027   , beam_origin(beamOrigin)
0028   , phi_diff_cut(phiDiffCut)
0029   , DCA_cut(DCACut)
0030   , N_clu_cutl(NCluCutl)
0031   , N_clu_cut(NCluCut)
0032   , angle_diff_new_l(angleDiffNew_l)
0033   , angle_diff_new_r(angleDiffNew_r)
0034   , peek(peekCut)
0035   , print_message_opt(printMessageOpt)
0036 {
0037   gErrorIgnoreLevel = kWarning;  // note : To not print the "print plot info."
0038 
0039   // Init();
0040   plot_text = (run_type == "MC") ? "Simulation" : "Work-in-progress";
0041 
0042   cluster_pair_vec.clear();
0043 
0044   current_vtxX = beam_origin.first;
0045   current_vtxY = beam_origin.second;
0046 }
0047 
0048 INTTXYvtx::~INTTXYvtx()
0049 {
0050   // InitHist
0051   for (auto& itr : m_v_hist)
0052   {
0053     //--std::cout<<"del : "<<itr->GetTitle()<<" "<<std::hex<<(long)itr<<std::hex<<std::endl;
0054     delete itr;
0055   }
0056   //--for(auto& itr: m_v_hist){
0057   //--  std::cout<<"after del : "<<std::hex<<(long)itr<<std::hex<<std::endl;
0058   //--}
0059 
0060   // InitGraph
0061   delete angle_diff_inner_phi_peak_profile_graph;
0062   delete angle_diff_outer_phi_peak_profile_graph;
0063   delete DCA_distance_inner_phi_peak_profile_graph;
0064   delete DCA_distance_outer_phi_peak_profile_graph;
0065 
0066   // created in FillLine_FindVertex
0067   delete xy_hist;
0068   delete xy_hist_bkgrm;
0069 
0070   // InitRest
0071   delete cos_fit;
0072   delete gaus_fit;
0073   delete horizontal_fit_inner;
0074   delete horizontal_fit_angle_diff_inner;
0075   delete horizontal_fit_outer;
0076   delete horizontal_fit_angle_diff_outer;
0077 
0078   delete c1;
0079   delete ltx;
0080   delete draw_text;
0081 }
0082 
0083 void INTTXYvtx::Init()
0084 {
0085   if (m_enable_drawhist || m_savehist)
0086   {
0087     if (!std::filesystem::exists(out_folder_directory))
0088     {
0089       std::filesystem::create_directory(out_folder_directory);
0090     }
0091   }
0092 
0093   InitHist();
0094   InitGraph();
0095   InitRest();
0096 
0097   //--    InitTreeOut();
0098 
0099   m_initialized = true;
0100 }
0101 
0102 void INTTXYvtx::InitGraph()
0103 {
0104   angle_diff_inner_phi_peak_profile_graph = new TGraph();
0105   angle_diff_outer_phi_peak_profile_graph = new TGraph();
0106   DCA_distance_inner_phi_peak_profile_graph = new TGraph();
0107   DCA_distance_outer_phi_peak_profile_graph = new TGraph();
0108 }
0109 
0110 void INTTXYvtx::InitRest()
0111 {
0112   horizontal_fit_inner = new TF1("horizontal_fit_inner", "pol0", 0, 360);
0113   horizontal_fit_inner->SetLineWidth(2);
0114   horizontal_fit_inner->SetLineColor(2);
0115 
0116   horizontal_fit_angle_diff_inner = new TF1("horizontal_fit_angle_diff_inner", "pol0", 0, 360);
0117   horizontal_fit_angle_diff_inner->SetLineWidth(2);
0118   horizontal_fit_angle_diff_inner->SetLineColor(2);
0119 
0120   horizontal_fit_outer = new TF1("horizontal_fit_outer", "pol0", 0, 360);
0121   horizontal_fit_outer->SetLineWidth(2);
0122   horizontal_fit_outer->SetLineColor(2);
0123 
0124   horizontal_fit_angle_diff_outer = new TF1("horizontal_fit_angle_diff_outer", "pol0", 0, 360);
0125   horizontal_fit_angle_diff_outer->SetLineWidth(2);
0126   horizontal_fit_angle_diff_outer->SetLineColor(2);
0127 
0128   if (m_enable_qa)
0129   {
0130     cos_fit = new TF1("cos_fit", cos_func, 0, 360, 4);
0131     cos_fit->SetParNames("[A]", "[B]", "[C]", "[D]");
0132     cos_fit->SetLineColor(2);
0133 
0134     gaus_fit = new TF1("gaus_fit", InttVertexUtil::gaus_func, 0, 360, 4);
0135     gaus_fit->SetLineColor(4);
0136     gaus_fit->SetLineWidth(1);
0137     gaus_fit->SetParNames("size", "mean", "width", "offset");
0138     gaus_fit->SetNpx(1000);
0139   }
0140 
0141   if (m_enable_drawhist)
0142   {
0143     c1 = new TCanvas("", "", 950, 800);
0144     c1->cd();
0145 
0146     ltx = new TLatex();
0147     ltx->SetNDC();
0148     ltx->SetTextSize(0.045);
0149     ltx->SetTextAlign(31);
0150 
0151     draw_text = new TLatex();
0152     draw_text->SetNDC();
0153     draw_text->SetTextSize(0.03);
0154   }
0155 }
0156 
0157 void INTTXYvtx::InitHist()
0158 {
0159   ///////////////////////////////////////
0160   // QA histograms in process_evt
0161   if (m_enable_qa)
0162   {
0163     N_cluster_correlation = new TH2F("N_cluster_correlation", "N_cluster_correlation", 100, 0, 4000, 100, 0, 4000);
0164     N_cluster_correlation->SetStats(false);
0165     N_cluster_correlation->GetXaxis()->SetTitle("inner N_cluster");
0166     N_cluster_correlation->GetYaxis()->SetTitle("Outer N_cluster");
0167     N_cluster_correlation->GetXaxis()->SetNdivisions(505);
0168     m_v_hist.push_back(N_cluster_correlation);
0169 
0170     N_cluster_correlation_close = new TH2F("N_cluster_correlation_close", "N_cluster_correlation_close", 100, 0, 500, 100, 0, 500);
0171     N_cluster_correlation_close->SetStats(false);
0172     N_cluster_correlation_close->GetXaxis()->SetTitle("inner N_cluster");
0173     N_cluster_correlation_close->GetYaxis()->SetTitle("Outer N_cluster");
0174     N_cluster_correlation_close->GetXaxis()->SetNdivisions(505);
0175     m_v_hist.push_back(N_cluster_correlation_close);
0176 
0177     inner_pos_xy = new TH2F("inner_pos_xy", "inner_pos_xy", 360, -100, 100, 360, -100, 100);
0178     inner_pos_xy->SetStats(false);
0179     inner_pos_xy->GetXaxis()->SetTitle("X axis [mm]");
0180     inner_pos_xy->GetYaxis()->SetTitle("Y axis [mm]");
0181     inner_pos_xy->GetXaxis()->SetNdivisions(505);
0182     m_v_hist.push_back(inner_pos_xy);
0183 
0184     outer_pos_xy = new TH2F("outer_pos_xy", "outer_pos_xy", 360, -150, 150, 360, -150, 150);
0185     outer_pos_xy->SetStats(false);
0186     outer_pos_xy->GetXaxis()->SetTitle("X axis [mm]");
0187     outer_pos_xy->GetYaxis()->SetTitle("Y axis [mm]");
0188     outer_pos_xy->GetXaxis()->SetNdivisions(505);
0189     m_v_hist.push_back(outer_pos_xy);
0190 
0191     // inner_outer_pos_xy = new TH2F("inner_outer_pos_xy","inner_outer_pos_xy",360,-150,150,360,-150,150);
0192     inner_outer_pos_xy = new TH2F("inner_outer_pos_xy", "inner_outer_pos_xy", 360, -150, 150, 360, -150, 150);
0193     inner_outer_pos_xy->SetStats(false);
0194     inner_outer_pos_xy->GetXaxis()->SetTitle("X axis [mm]");
0195     inner_outer_pos_xy->GetYaxis()->SetTitle("Y axis [mm]");
0196     inner_outer_pos_xy->GetXaxis()->SetNdivisions(505);
0197     m_v_hist.push_back(inner_outer_pos_xy);
0198   }
0199 
0200   ///////////////////////////////////////
0201   // histograms for quadrant method
0202   DCA_distance_inner_phi = new TH2F("DCA_distance_inner_phi", "DCA_distance_inner_phi", 100, 0, 360, 100, -10, 10);
0203   DCA_distance_inner_phi->SetStats(false);
0204   DCA_distance_inner_phi->GetXaxis()->SetTitle("Inner phi [degree]");
0205   DCA_distance_inner_phi->GetYaxis()->SetTitle("DCA [mm]");
0206   DCA_distance_inner_phi->GetXaxis()->SetNdivisions(505);
0207   m_v_hist.push_back(DCA_distance_inner_phi);
0208 
0209   angle_diff_inner_phi = new TH2F("angle_diff_inner_phi", "angle_diff_inner_phi", 361, 0, 361, 100, -1.5, 1.5);
0210   angle_diff_inner_phi->SetStats(false);
0211   angle_diff_inner_phi->GetXaxis()->SetTitle("Inner phi [degree]");
0212   angle_diff_inner_phi->GetYaxis()->SetTitle("Inner - Outer [degree]");
0213   angle_diff_inner_phi->GetXaxis()->SetNdivisions(505);
0214   m_v_hist.push_back(angle_diff_inner_phi);
0215 
0216   ///////////////////////////////////////
0217   // note : it's for the geometry correction
0218   angle_diff_inner_phi_peak_final = new TH2F("angle_diff_inner_phi_peak_final",
0219                                              "angle_diff_inner_phi_peak_final", 361, 0, 361, 100, -1.5, 1.5);
0220   angle_diff_inner_phi_peak_final->SetStats(false);
0221   angle_diff_inner_phi_peak_final->GetXaxis()->SetTitle("Inner phi [degree]");
0222   angle_diff_inner_phi_peak_final->GetYaxis()->SetTitle("Inner - Outer [degree]");
0223   angle_diff_inner_phi_peak_final->GetXaxis()->SetNdivisions(505);
0224   m_v_hist.push_back(angle_diff_inner_phi_peak_final);
0225 
0226   DCA_distance_inner_phi_peak_final = new TH2F("DCA_distance_inner_phi_peak_final",
0227                                                "DCA_distance_inner_phi_peak_final", 100, 0, 360, 100, -10, 10);
0228   DCA_distance_inner_phi_peak_final->SetStats(false);
0229   DCA_distance_inner_phi_peak_final->GetXaxis()->SetTitle("Inner phi [degree]");
0230   DCA_distance_inner_phi_peak_final->GetYaxis()->SetTitle("DCA [mm]");
0231   DCA_distance_inner_phi_peak_final->GetXaxis()->SetNdivisions(505);
0232   m_v_hist.push_back(DCA_distance_inner_phi_peak_final);
0233 
0234   // QA
0235   angle_diff = new TH1F("angle_diff", "angle_diff", 100, 0, 5);
0236   angle_diff->SetStats(false);
0237   angle_diff->GetXaxis()->SetTitle("|Inner - Outer| [degree]");
0238   angle_diff->GetYaxis()->SetTitle("Entry");
0239   angle_diff->GetXaxis()->SetNdivisions(505);
0240   m_v_hist.push_back(angle_diff);
0241 
0242   angle_diff_new = new TH1F("angle_diff_new", "angle_diff_new", 100, angle_diff_new_l, angle_diff_new_r);
0243   angle_diff_new->SetStats(false);
0244   angle_diff_new->GetXaxis()->SetTitle("|Inner - Outer| [degree]");
0245   angle_diff_new->GetYaxis()->SetTitle("Entry");
0246   angle_diff_new->GetXaxis()->SetNdivisions(505);
0247   m_v_hist.push_back(angle_diff_new);
0248 
0249   DCA_distance = new TH1F("DCA_distance", "DCA_distance", 100, -10, 10);
0250   DCA_distance->SetStats(false);
0251   DCA_distance->GetXaxis()->SetTitle("DCA [mm]");
0252   DCA_distance->GetYaxis()->SetTitle("Entry");
0253   DCA_distance->GetXaxis()->SetNdivisions(505);
0254   m_v_hist.push_back(DCA_distance);
0255 
0256   if (m_enable_qa)
0257   {
0258     DCA_distance_outer_phi = new TH2F("DCA_distance_outer_phi", "DCA_distance_outer_phi", 100, 0, 360, 100, -10, 10);
0259     DCA_distance_outer_phi->SetStats(false);
0260     DCA_distance_outer_phi->GetXaxis()->SetTitle("Outer phi [degree]");
0261     DCA_distance_outer_phi->GetYaxis()->SetTitle("DCA [mm]");
0262     DCA_distance_outer_phi->GetXaxis()->SetNdivisions(505);
0263     m_v_hist.push_back(DCA_distance_outer_phi);
0264 
0265     angle_diff_outer_phi = new TH2F("angle_diff_outer_phi", "angle_diff_outer_phi", 361, 0, 361, 100, -1.5, 1.5);
0266     angle_diff_outer_phi->SetStats(false);
0267     angle_diff_outer_phi->GetXaxis()->SetTitle("Outer phi [degree]");
0268     angle_diff_outer_phi->GetYaxis()->SetTitle("Inner - Outer [degree]");
0269     angle_diff_outer_phi->GetXaxis()->SetNdivisions(505);
0270     m_v_hist.push_back(angle_diff_outer_phi);
0271 
0272     angle_correlation = new TH2F("angle_correlation", "angle_correlation", 361, 0, 361, 361, 0, 361);
0273     angle_correlation->SetStats(false);
0274     angle_correlation->GetXaxis()->SetTitle("Inner Phi [degree]");
0275     angle_correlation->GetYaxis()->SetTitle("Outer Phi [degree]");
0276     angle_correlation->GetXaxis()->SetNdivisions(505);
0277     m_v_hist.push_back(angle_correlation);
0278 
0279     angle_diff_DCA_dist = new TH2F("angle_diff_DCA_dist", "angle_diff_DCA_dist", 100, -1.5, 1.5, 100, -3.5, 3.5);
0280     angle_diff_DCA_dist->SetStats(false);
0281     angle_diff_DCA_dist->GetXaxis()->SetTitle("Inner - Outer [degree]");
0282     angle_diff_DCA_dist->GetYaxis()->SetTitle("DCA distance [mm]");
0283     angle_diff_DCA_dist->GetXaxis()->SetNdivisions(505);
0284     m_v_hist.push_back(angle_diff_DCA_dist);
0285 
0286     DCA_point = new TH2F("DCA_point", "DCA_point", 100, -10, 10, 100, -10, 10);
0287     DCA_point->SetStats(false);
0288     DCA_point->GetXaxis()->SetTitle("X pos [mm]");
0289     DCA_point->GetYaxis()->SetTitle("Y pos [mm]");
0290     DCA_point->GetXaxis()->SetNdivisions(505);
0291     m_v_hist.push_back(DCA_point);
0292 
0293     DCA_distance_inner_X = new TH2F("DCA_distance_inner_X", "DCA_distance_inner_X", 100, -100, 100, 100, -10, 10);
0294     DCA_distance_inner_X->SetStats(false);
0295     DCA_distance_inner_X->GetXaxis()->SetTitle("Inner cluster X [mm]");
0296     DCA_distance_inner_X->GetYaxis()->SetTitle("DCA [mm]");
0297     DCA_distance_inner_X->GetXaxis()->SetNdivisions(505);
0298     m_v_hist.push_back(DCA_distance_inner_X);
0299 
0300     DCA_distance_inner_Y = new TH2F("DCA_distance_inner_Y", "DCA_distance_inner_Y", 100, -100, 100, 100, -10, 10);
0301     DCA_distance_inner_Y->SetStats(false);
0302     DCA_distance_inner_Y->GetXaxis()->SetTitle("Inner cluster Y [mm]");
0303     DCA_distance_inner_Y->GetYaxis()->SetTitle("DCA [mm]");
0304     DCA_distance_inner_Y->GetXaxis()->SetNdivisions(505);
0305     m_v_hist.push_back(DCA_distance_inner_Y);
0306 
0307     DCA_distance_outer_X = new TH2F("DCA_distance_outer_X", "DCA_distance_outer_X", 100, -130, 130, 100, -10, 10);
0308     DCA_distance_outer_X->SetStats(false);
0309     DCA_distance_outer_X->GetXaxis()->SetTitle("Outer cluster X [mm]");
0310     DCA_distance_outer_X->GetYaxis()->SetTitle("DCA [mm]");
0311     DCA_distance_outer_X->GetXaxis()->SetNdivisions(505);
0312     m_v_hist.push_back(DCA_distance_outer_X);
0313 
0314     DCA_distance_outer_Y = new TH2F("DCA_distance_outer_Y", "DCA_distance_outer_Y", 100, -130, 130, 100, -10, 10);
0315     DCA_distance_outer_Y->SetStats(false);
0316     DCA_distance_outer_Y->GetXaxis()->SetTitle("Outer cluster Y [mm]");
0317     DCA_distance_outer_Y->GetYaxis()->SetTitle("DCA [mm]");
0318     DCA_distance_outer_Y->GetXaxis()->SetNdivisions(505);
0319     m_v_hist.push_back(DCA_distance_outer_Y);
0320 
0321     ///////////////////////////////////////
0322     // note : it's for the geometry correction
0323     angle_diff_outer_phi_peak_final = new TH2F("angle_diff_outer_phi_peak_final",
0324                                                "angle_diff_outer_phi_peak_final", 361, 0, 361, 100, -1.5, 1.5);
0325     angle_diff_outer_phi_peak_final->SetStats(false);
0326     angle_diff_outer_phi_peak_final->GetXaxis()->SetTitle("Outer phi [degree]");
0327     angle_diff_outer_phi_peak_final->GetYaxis()->SetTitle("Inner - Outer [degree]");
0328     angle_diff_outer_phi_peak_final->GetXaxis()->SetNdivisions(505);
0329     m_v_hist.push_back(angle_diff_outer_phi_peak_final);
0330 
0331     DCA_distance_outer_phi_peak_final = new TH2F("DCA_distance_outer_phi_peak_final",
0332                                                  "DCA_distance_outer_phi_peak_final", 100, 0, 360, 100, -10, 10);
0333     DCA_distance_outer_phi_peak_final->SetStats(false);
0334     DCA_distance_outer_phi_peak_final->GetXaxis()->SetTitle("Outer phi [degree]");
0335     DCA_distance_outer_phi_peak_final->GetYaxis()->SetTitle("DCA [mm]");
0336     DCA_distance_outer_phi_peak_final->GetXaxis()->SetNdivisions(505);
0337     m_v_hist.push_back(DCA_distance_outer_phi_peak_final);
0338 
0339     angle_diff_new_bkg_remove_final = new TH1F("angle_diff_new_bkg_remove_final",
0340                                                "angle_diff_new_bkg_remove_final", 100, angle_diff_new_l, angle_diff_new_r);
0341     angle_diff_new_bkg_remove_final->SetStats(false);
0342     angle_diff_new_bkg_remove_final->GetXaxis()->SetTitle("|Inner - Outer| [degree]");
0343     angle_diff_new_bkg_remove_final->GetYaxis()->SetTitle("Entry");
0344     angle_diff_new_bkg_remove_final->GetXaxis()->SetNdivisions(505);
0345     m_v_hist.push_back(angle_diff_new_bkg_remove_final);
0346   }
0347 }
0348 
0349 // note : this function only prepare the pairs for the vertex XY calculation, it's like a general vertex for the whole run
0350 void INTTXYvtx::ProcessEvt(
0351     int event_i,
0352     const std::vector<clu_info>& temp_sPH_inner_nocolumn_vec,
0353     const std::vector<clu_info>& temp_sPH_outer_nocolumn_vec,
0354     const std::vector<std::vector<double>>& /*temp_sPH_nocolumn_vec*/,
0355     const std::vector<std::vector<double>>& /*temp_sPH_nocolumn_rz_vec*/,
0356     int NvtxMC,
0357     double /*TrigZvtxMC*/,
0358     bool PhiCheckTag,
0359     Long64_t /*bco_full*/
0360 )
0361 {
0362   if (!m_initialized)
0363   {
0364     std::cout << "INTTXYvtx is not initialized, abort in ProcessEvt" << std::endl;
0365     exit(1);
0366   }
0367 
0368   if (print_message_opt && event_i % 10000 == 0)
0369   {
0370     std::cout << "In INTTXYvtx class, running event : " << event_i << std::endl;
0371   }
0372 
0373   total_NClus = temp_sPH_inner_nocolumn_vec.size() + temp_sPH_outer_nocolumn_vec.size();
0374 
0375   // note : the Move these two in the beginning of the function, in order to avoid those event-reject cuts
0376   if (m_enable_qa)
0377   {
0378     N_cluster_correlation->Fill(temp_sPH_inner_nocolumn_vec.size(), temp_sPH_outer_nocolumn_vec.size());
0379     N_cluster_correlation_close->Fill(temp_sPH_inner_nocolumn_vec.size(), temp_sPH_outer_nocolumn_vec.size());
0380   }
0381 
0382   if (total_NClus < zvtx_cal_require)
0383   {
0384     return;
0385     std::cout << "return confirmation" << std::endl;
0386   }
0387 
0388   if (run_type == "MC" && NvtxMC != 1)
0389   {
0390     return;
0391     std::cout << "In INTTXYvtx class, event : " << event_i
0392               << " Nvtx : " << NvtxMC << " Nvtx more than one " << std::endl;
0393   }
0394   if (PhiCheckTag == false)
0395   {
0396     return;
0397     std::cout << "In INTTXYvtx class, event : " << event_i
0398               << " Nvtx : " << NvtxMC << " Not full phi has hits " << std::endl;
0399   }
0400 
0401   if (temp_sPH_inner_nocolumn_vec.size() < 10 || temp_sPH_outer_nocolumn_vec.size() < 10 || total_NClus > N_clu_cut || total_NClus < N_clu_cutl)
0402   {
0403     return;
0404   }
0405 
0406   //-------------------------------
0407   // tracklet reconstruction accumulated multiple events
0408   for (auto& inner_i : temp_sPH_inner_nocolumn_vec)
0409   {
0410     for (auto& outer_i : temp_sPH_outer_nocolumn_vec)
0411     {
0412       // note : try to ease the analysis and also make it quick.
0413       if (fabs(inner_i.phi - outer_i.phi) < 7)  // todo : the pre phi cut is here, can be optimized
0414       {
0415         cluster_pair_vec.push_back({{inner_i.x,
0416                                      inner_i.y},
0417                                     {outer_i.x,
0418                                      outer_i.y}});
0419       }
0420     }
0421   }
0422 
0423   //-------------------------------
0424   // QA histogram
0425   if (m_enable_qa)
0426   {
0427     for (auto& inner_i : temp_sPH_inner_nocolumn_vec)
0428     {
0429       inner_pos_xy->Fill(inner_i.x, inner_i.y);
0430       inner_outer_pos_xy->Fill(inner_i.x, inner_i.y);
0431     }
0432 
0433     for (auto& outer_i : temp_sPH_outer_nocolumn_vec)
0434     {
0435       outer_pos_xy->Fill(outer_i.x, outer_i.y);
0436       inner_outer_pos_xy->Fill(outer_i.x, outer_i.y);
0437     }
0438   }
0439 
0440   //--  std::cout<<"  "<<event_i<<" clusterpair:size : "<<cluster_pair_vec.size()<<std::endl;
0441 }
0442 
0443 void INTTXYvtx::ClearEvt()
0444 {
0445   return;
0446 }
0447 
0448 unsigned long INTTXYvtx::GetVecNele()
0449 {
0450   return cluster_pair_vec.size();
0451 }
0452 
0453 std::pair<double, double> INTTXYvtx::GetFinalVTXxy()
0454 {
0455   return {current_vtxX, current_vtxY};
0456 }
0457 
0458 std::pair<std::vector<TH2F*>, std::vector<TH1F*>> INTTXYvtx::GetHistFinal()
0459 {
0460   return {
0461       {DCA_distance_inner_phi_peak_final,
0462        angle_diff_inner_phi_peak_final,
0463        DCA_distance_outer_phi_peak_final,
0464        angle_diff_outer_phi_peak_final,
0465        xy_hist,
0466        xy_hist_bkgrm},
0467       {angle_diff_new_bkg_remove_final}};
0468 }
0469 
0470 void INTTXYvtx::PrintPlots()
0471 {
0472   if (!m_initialized)
0473   {
0474     std::cout << "INTTXYvtx is not initialized, abort in PrintPlots" << std::endl;
0475     exit(1);
0476   }
0477 
0478   if (m_enable_drawhist && m_enable_qa)
0479   {
0480     std::string s_inttlabel = (boost::format("#it{#bf{sPHENIX INTT}} %s") % plot_text).str().c_str();
0481     // note : -----------------------------------------------------------------------------------------
0482     inner_outer_pos_xy->Draw("colz0");
0483     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
0484     c1->Print((boost::format("%s/xyvtx_qa.pdf(") % out_folder_directory).str().c_str());
0485     c1->Clear();
0486 
0487     // note : -----------------------------------------------------------------------------------------
0488     inner_pos_xy->Draw("colz0");
0489     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
0490     c1->Print((boost::format("%s/xyvtx_qa.pdf") % out_folder_directory).str().c_str());
0491     c1->Clear();
0492 
0493     // note : -----------------------------------------------------------------------------------------
0494     outer_pos_xy->Draw("colz0");
0495     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
0496     c1->Print((boost::format("%s/xyvtx_qa.pdf") % out_folder_directory).str().c_str());
0497     c1->Clear();
0498 
0499     // note : -----------------------------------------------------------------------------------------
0500     N_cluster_correlation->Draw("colz0");
0501     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
0502     c1->Print((boost::format("%s/xyvtx_qa.pdf") % out_folder_directory).str().c_str());
0503     c1->Clear();
0504 
0505     // note : -----------------------------------------------------------------------------------------
0506     N_cluster_correlation_close->Draw("colz0");
0507     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
0508     c1->Print((boost::format("%s/xyvtx_qa.pdf)") % out_folder_directory).str().c_str());
0509     c1->Clear();
0510   }
0511 }
0512 
0513 //--------------------------------------------------------------------------
0514 // quadrant method
0515 std::vector<std::pair<double, double>> INTTXYvtx::MacroVTXSquare(double length, int N_trial)
0516 {
0517   if (!m_initialized)
0518   {
0519     std::cout << "INTTXYvtx is not initialized, abort in MacroVTXSquare" << std::endl;
0520     exit(1);
0521   }
0522 
0523   bool draw_plot_opt = m_enable_drawhist;
0524 
0525   const double original_length = length;
0526   std::pair<double, double> origin = {0, 0};
0527   std::vector<std::pair<double, double>> vtx_vec = Get4vtx(origin, length);  // vtx_vec.push_back(origin);
0528 
0529   int small_index{0};
0530   std::vector<double> small_info_vec(18, -999);
0531   std::vector<double> grr_x{};
0532   std::vector<double> grr_E{};
0533   std::vector<double> grr_y{};
0534 
0535   std::vector<double> All_FitError_DCA_Y{};
0536   std::vector<double> All_FitError_DCA_X{};
0537   std::vector<double> All_FitError_angle_Y{};
0538   std::vector<double> All_FitError_angle_X{};
0539 
0540   std::vector<double> Winner_FitError_DCA_Y{};
0541   std::vector<double> Winner_FitError_DCA_X{};
0542   std::vector<double> Winner_FitError_angle_Y{};
0543   std::vector<double> Winner_FitError_angle_X{};
0544 
0545   if (print_message_opt == true)
0546   {
0547     std::cout << "In INTTXYvtx::MacroVTXSquare, N pairs : " << cluster_pair_vec.size() << std::endl;
0548 
0549     std::cout << N_trial << " runs, smart. which gives you the resolution down to "
0550               << length / pow(2, N_trial) << " mm" << std::endl;
0551   }
0552 
0553   if (cluster_pair_vec.size() == 0)
0554   {  // minimum tracklet cut. need to be tuned
0555     return {
0556         beam_origin,  // note : the best vertex
0557         origin,       // note : the origin in that trial
0558         {0, 0},       // note : horizontal_fit_inner -> GetParError(0),  horizontal_fit_angle_diff_inner -> GetParError(0)
0559         {0, 0},       // note : horizontal_fit_inner -> GetParameter(0), horizontal_fit_angle_diff_inner -> GetParameter(0)
0560         {0, 0},       // note : horizontal_fit_outer -> GetParError(0),  horizontal_fit_angle_diff_outer -> GetParError(0)
0561         {0, 0},       // note : horizontal_fit_outer -> GetParameter(0), horizontal_fit_angle_diff_outer -> GetParameter(0)
0562         {0, 0},       // note : the mean and stddev of angle_diff
0563         {0, 0},       // note : the mean and stddev of DCA_distance
0564         {0, 0},       // note : the mean and stddev of angle_diff, but with the background removed
0565     };
0566   }
0567 
0568   //-------------------------------
0569   // info_vec contents
0570   //   0 : angle_diff_new_bkg_remove       -> Mean,
0571   //   1 : angle_diff_new_bkg_remove       -> StdDev,      // note : angle diff stddev and error (1D histogram)
0572   //   2 : horizontal_fit_inner            -> Chisq / NDF),
0573   //   3 : horizontal_fit_inner            -> ParErr(0),   // note : inner DCA, pol0
0574   //   4 : horizontal_fit_angle_diff_inner -> Chisq / NDF),
0575   //   5 : horizontal_fit_angle_diff_inner -> ParErr(0),   // note : inner angle diff, pol0
0576   //   6 : horizontal_fit_outer            -> Chisq / NDF),
0577   //   7 : horizontal_fit_outer            -> ParErr(0),   // note : outer DCA, pol0
0578   //   8 : horizontal_fit_angle_diff_outer -> Chisq / NDF),
0579   //   9 : horizontal_fit_angle_diff_outer -> ParErr(0),   // note : outer angle diff, pol0
0580   //  10 : horizontal_fit_inner            -> Par(0),
0581   //  11 : horizontal_fit_angle_diff_inner -> Par(0),
0582   //  12 : horizontal_fit_outer            -> Par(0),
0583   //  13 : horizontal_fit_angle_diff_outer -> Par(0),
0584   //  14 : angle_diff                      -> Mean,
0585   //  15 : angle_diff                      -> StdDev,
0586   //  16 : DCA_distance                    -> Mean,
0587   //  17 : DCA_distance                    -> StdDev
0588   //
0589 
0590   if (m_enable_drawhist)
0591   {
0592     c1->cd();
0593     c1->Range(0, 0, 1, 1);
0594     ltx->DrawLatex(0.5, 0.5, "QA plots for quadrant method");
0595     c1->Print((boost::format("%s/%s(") % out_folder_directory % m_quad_pdfname).str().c_str());
0596     c1->Clear();
0597   }
0598 
0599   // current algorithm uses info_vec[3] and [5], others are for QA
0600   for (int i = 0; i < N_trial; i++)
0601   {
0602     if (print_message_opt == true)
0603     {
0604       std::cout << "~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~"
0605                 << " step " << i << " "
0606                 << "~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~" << std::endl;
0607     }
0608     for (unsigned int i1 = 0; i1 < vtx_vec.size(); i1++)
0609     {
0610       if (print_message_opt == true)
0611       {
0612         std::cout << "tested vertex : " << vtx_vec[i1].first << " " << vtx_vec[i1].second << std::endl;
0613       }
0614 
0615       if (m_enable_drawhist)
0616       {
0617         c1->cd();
0618         c1->Range(0, 0, 1, 1);
0619         ltx->DrawLatex(0.5, 0.5, (boost::format("New_trial_square_%i_%i") % i % i1).str().c_str());
0620         c1->Print((boost::format("%s/%s") % out_folder_directory.c_str() % m_quad_pdfname).str().c_str());
0621         c1->Clear();
0622       }
0623 
0624       current_vtxX = vtx_vec[i1].first;
0625       current_vtxY = vtx_vec[i1].second;
0626 
0627       std::vector<double> info_vec = subMacroVTXxyCorrection(i, i1, draw_plot_opt);
0628 
0629       if (print_message_opt == true)
0630       {
0631         std::cout << "trial : " << i
0632                   << " vertex : " << i1
0633                   << " DCA fit error : " << info_vec[3]
0634                   << " angle diff fit error : " << info_vec[5] << std::endl;
0635       }
0636 
0637       All_FitError_DCA_Y.push_back(info_vec[3]);
0638       All_FitError_DCA_X.push_back(i);
0639       All_FitError_angle_Y.push_back(info_vec[5]);
0640       All_FitError_angle_X.push_back(i);
0641 
0642       if ((i1 == 0) ||
0643           (info_vec[3] < small_info_vec[3] &&
0644            info_vec[5] < small_info_vec[5])  // note : the fit error of the pol0 fit
0645       )
0646       {
0647         small_info_vec = info_vec;
0648         small_index = i1;
0649 
0650         TH2F_FakeClone(DCA_distance_inner_phi_peak, DCA_distance_inner_phi_peak_final);
0651         TH2F_FakeClone(angle_diff_inner_phi_peak, angle_diff_inner_phi_peak_final);
0652         if (m_enable_qa)
0653         {
0654           TH2F_FakeClone(DCA_distance_outer_phi_peak, DCA_distance_outer_phi_peak_final);
0655           TH2F_FakeClone(angle_diff_outer_phi_peak, angle_diff_outer_phi_peak_final);
0656           TH1F_FakeClone(angle_diff_new_bkg_remove, angle_diff_new_bkg_remove_final);
0657         }
0658       }
0659       if (print_message_opt == true)
0660       {
0661         std::cout << " " << std::endl;
0662       }
0663 
0664       ClearHist(1);
0665     }
0666 
0667     if (print_message_opt == true)
0668     {
0669       std::cout << "the Quadrant " << small_index << " won the competition" << std::endl;
0670     }
0671 
0672     Winner_FitError_DCA_Y.push_back(small_info_vec[3]);
0673     Winner_FitError_DCA_X.push_back(i);
0674     Winner_FitError_angle_Y.push_back(small_info_vec[5]);
0675     Winner_FitError_angle_X.push_back(i);
0676 
0677     grr_x.push_back(i);
0678     grr_y.push_back(small_index);
0679     grr_E.push_back(0);
0680 
0681     // note : generating the new 4 vertex for the next comparison
0682     // note : start to shrink the square
0683     if (i != N_trial - 1)
0684     {
0685       origin = {(vtx_vec[small_index].first + origin.first) / 2.,
0686                 (vtx_vec[small_index].second + origin.second) / 2.};
0687 
0688       // std::cout<<"test : "<<origin.first<<" "<<origin.second<<" length: "<<length<<std::endl;
0689       // if (small_index == 4) {length /= 1.5;}
0690       // else {length /= 2.;}
0691       length /= 2.;
0692       vtx_vec = Get4vtx(origin, length);  // vtx_vec.push_back(origin);
0693     }
0694   }
0695 
0696   if (m_enable_drawhist)
0697   {
0698     c1->cd();
0699     c1->Print((boost::format("%s/%s)") % out_folder_directory.c_str() % m_quad_pdfname).str().c_str());
0700     c1->Clear();
0701   }
0702 
0703   if (draw_plot_opt == true)
0704   {
0705     DrawTGraphErrors(grr_x, grr_y, grr_E, grr_E, out_folder_directory,
0706                      {
0707                          (boost::format("Square_scan_history_%.1fmm_%iTrials") % original_length % N_trial).str().c_str()  // title
0708                          ,
0709                          "nth scan"  // x_title
0710                          ,
0711                          "Winner index"  // y_title
0712                          ,
0713                          "APL"  // draw option
0714                          ,
0715                          "quadorant_qa.pdf("  // pdf name
0716                      });
0717     Draw2TGraph(All_FitError_angle_X, All_FitError_angle_Y, Winner_FitError_angle_X, Winner_FitError_angle_Y, out_folder_directory,
0718                 {
0719                     (boost::format("Angle_diff_fit_error_%iTrials") % N_trial).str().c_str()  // title
0720                     ,
0721                     "n iteration"  // x_title
0722                     ,
0723                     "#Delta#phi fit error [degree]"  // y_title
0724                     ,
0725                     ""  // draw option    // draw option
0726                     ,
0727                     "quadorant_qa.pdf"  // pdf name   // pdf name
0728                 });
0729     Draw2TGraph(All_FitError_DCA_X, All_FitError_DCA_Y, Winner_FitError_DCA_X, Winner_FitError_DCA_Y, out_folder_directory,
0730                 {
0731                     (boost::format("DCA_fit_error_%iTrials") % N_trial).str().c_str()  // title
0732                     ,
0733                     "n iteration"  // x_title
0734                     ,
0735                     "DCA fit error [mm]"  // y_title
0736                     ,
0737                     ""  // draw option
0738                     ,
0739                     "quadorant_qa.pdf)"  // pdf name
0740                 });
0741   }
0742 
0743   return {
0744       vtx_vec[small_index],                      // note : the best vertex
0745       origin,                                    // note : the origin in that trial
0746       {small_info_vec[3], small_info_vec[5]},    // note : horizontal_fit_inner -> GetParError(0),  horizontal_fit_angle_diff_inner -> GetParError(0)
0747       {small_info_vec[10], small_info_vec[11]},  // note : horizontal_fit_inner -> GetParameter(0), horizontal_fit_angle_diff_inner -> GetParameter(0)
0748       {small_info_vec[7], small_info_vec[9]},    // note : horizontal_fit_outer -> GetParError(0),  horizontal_fit_angle_diff_outer -> GetParError(0)
0749       {small_info_vec[12], small_info_vec[13]},  // note : horizontal_fit_outer -> GetParameter(0), horizontal_fit_angle_diff_outer -> GetParameter(0)
0750       {small_info_vec[14], small_info_vec[15]},  // note : the mean and stddev of angle_diff
0751       {small_info_vec[16], small_info_vec[17]},  // note : the mean and stddev of DCA_distance
0752       {small_info_vec[0], small_info_vec[1]},    // note : the mean and stddev of angle_diff, but with the background removed
0753   };
0754 }
0755 
0756 std::vector<double> INTTXYvtx::subMacroVTXxyCorrection(int test_index, int trial_index, bool draw_plot_opt)
0757 {
0758   int true_trial_index = test_index * 4 + trial_index;
0759   std::vector<double> out_vec = GetVTXxyCorrection_new(true_trial_index);
0760 
0761   std::string sub_out_folder_name{};
0762   if (draw_plot_opt == true)
0763   {
0764     sub_out_folder_name = (boost::format("%s/New_trial_square_%i_%i") % out_folder_directory % test_index % trial_index).str();
0765 
0766     //--if (std::filesystem::exists(sub_out_folder_name.c_str()) == false) {
0767     //--  system((boost::format("mkdir %s",sub_out_folder_name.c_str()));
0768     //--}
0769 
0770     // PrintPlotsVTXxy(sub_out_folder_name);
0771     PrintPlotsVTXxy();
0772   }
0773 
0774   return out_vec;
0775 }
0776 
0777 // note : {circle radius, possible correction angle, the chi2/NDF of pol0 fit}
0778 std::vector<double> INTTXYvtx::GetVTXxyCorrection_new(int trial_index)
0779 {
0780   if (print_message_opt == true)
0781   {
0782     std::cout << "Trial : " << trial_index
0783               << "---------------------------- ---------------------------- ----------------------------" << std::endl;
0784     std::cout << "Given vertex: " << current_vtxX << " " << current_vtxY << std::endl;
0785   }
0786 
0787   if (m_enable_qa)
0788   {
0789     cos_fit->SetParameters(4, 1.49143e-02, 185, 0.3);  // todo : we may have to apply more constaints on the fitting
0790     cos_fit->SetParLimits(2, 0, 360);                  // note : the peak location has to be positive
0791 
0792     // note : here is the test with a gaus fitting to find the peak
0793     gaus_fit->SetParameters(-4.5, 197, 50, 0);
0794     gaus_fit->SetParLimits(0, -100, 0);  // note : the gaus distribution points down
0795     // DCA_distance_inner_phi_peak_profile -> Fit(gaus_fit, "N","",100, 260);
0796     // std::cout<<"test, gaus fit range : "<<gaus_fit->GetParameter(1) - 25<<" "<<gaus_fit->GetParameter(1) + 25<<std::endl;
0797   }
0798 
0799   subMacroPlotWorking(true, 100, 260, 25);
0800 
0801   return {
0802       angle_diff_new_bkg_remove->GetMean(),
0803       angle_diff_new_bkg_remove->GetStdDev(),  // note : angle diff stddev and error (1D histogram)
0804       horizontal_fit_inner->GetChisquare() / double(horizontal_fit_inner->GetNDF()),
0805       horizontal_fit_inner->GetParError(0),  // note : inner DCA, pol0
0806       horizontal_fit_angle_diff_inner->GetChisquare() / double(horizontal_fit_angle_diff_inner->GetNDF()),
0807       horizontal_fit_angle_diff_inner->GetParError(0),  // note : inner angle diff, pol0
0808       horizontal_fit_outer->GetChisquare() / double(horizontal_fit_outer->GetNDF()),
0809       horizontal_fit_outer->GetParError(0),  // note : outer DCA, pol0
0810       horizontal_fit_angle_diff_outer->GetChisquare() / double(horizontal_fit_angle_diff_outer->GetNDF()),
0811       horizontal_fit_angle_diff_outer->GetParError(0),  // note : outer angle diff, pol0
0812       horizontal_fit_inner->GetParameter(0),
0813       horizontal_fit_angle_diff_inner->GetParameter(0),  // note : 10, 11
0814       horizontal_fit_outer->GetParameter(0),
0815       horizontal_fit_angle_diff_outer->GetParameter(0),  // note : 12, 13
0816       angle_diff->GetMean(),
0817       angle_diff->GetStdDev(),  // note : 14, 15
0818       DCA_distance->GetMean(),
0819       DCA_distance->GetStdDev(),  // note : 16, 17
0820   };
0821 }
0822 
0823 void INTTXYvtx::subMacroPlotWorking(
0824     bool phi_correction,
0825     double cos_fit_rangel,
0826     double cos_fit_ranger,
0827     double guas_fit_range)
0828 {
0829   //   3 : horizontal_fit_inner            -> ParErr(0),   // note : inner DCA, pol0
0830   //   5 : horizontal_fit_angle_diff_inner -> ParErr(0),   // note : inner angle diff, pol0
0831 
0832   for (auto& i : cluster_pair_vec)
0833   {
0834     std::vector<double> DCA_info_vec = calculateDistanceAndClosestPoint(
0835         i.first.x, i.first.y,
0836         i.second.x, i.second.y,
0837         current_vtxX, current_vtxY);
0838 
0839     double DCA_sign = calculateAngleBetweenVectors(
0840         i.second.x, i.second.y,
0841         i.first.x, i.first.y,
0842         current_vtxX, current_vtxY);
0843 
0844     if (phi_correction == true)
0845     {
0846       // std::cout<<"option selected "<<std::endl;
0847       Clus_InnerPhi_Offset = (i.first.y - current_vtxY < 0)
0848                                  ? atan2(i.first.y - current_vtxY, i.first.x - current_vtxX) * (180. / M_PI) + 360
0849                                  : atan2(i.first.y - current_vtxY, i.first.x - current_vtxX) * (180. / M_PI);
0850       Clus_OuterPhi_Offset = (i.second.y - current_vtxY < 0)
0851                                  ? atan2(i.second.y - current_vtxY, i.second.x - current_vtxX) * (180. / M_PI) + 360
0852                                  : atan2(i.second.y - current_vtxY, i.second.x - current_vtxX) * (180. / M_PI);
0853     }
0854     else  // note : phi_correction == false
0855     {
0856       Clus_InnerPhi_Offset = (i.first.y < 0)
0857                                  ? atan2(i.first.y, i.first.x) * (180. / M_PI) + 360
0858                                  : atan2(i.first.y, i.first.x) * (180. / M_PI);
0859       Clus_OuterPhi_Offset = (i.second.y < 0)
0860                                  ? atan2(i.second.y, i.second.x) * (180. / M_PI) + 360
0861                                  : atan2(i.second.y, i.second.x) * (180. / M_PI);
0862     }
0863 
0864     //----------------------
0865     // this is used for quadrant method
0866     DCA_distance_inner_phi->Fill(Clus_InnerPhi_Offset, DCA_sign);
0867     angle_diff_inner_phi->Fill(Clus_InnerPhi_Offset, Clus_InnerPhi_Offset - Clus_OuterPhi_Offset);
0868 
0869     //----------------------
0870     if (m_enable_qa)
0871     {
0872       DCA_distance_outer_phi->Fill(Clus_OuterPhi_Offset, DCA_sign);
0873       angle_diff_outer_phi->Fill(Clus_OuterPhi_Offset, Clus_InnerPhi_Offset - Clus_OuterPhi_Offset);
0874 
0875       angle_diff->Fill(std::abs(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset));
0876       angle_diff_new->Fill(std::abs(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset));
0877       DCA_distance->Fill(DCA_sign);
0878 
0879       // draw only
0880       angle_correlation->Fill(Clus_InnerPhi_Offset, Clus_OuterPhi_Offset);
0881       angle_diff_DCA_dist->Fill(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset, DCA_sign);
0882       DCA_point->Fill(DCA_info_vec[1], DCA_info_vec[2]);
0883 
0884       DCA_distance_inner_X->Fill(i.first.x, DCA_sign);
0885       DCA_distance_inner_Y->Fill(i.first.y, DCA_sign);
0886       DCA_distance_outer_X->Fill(i.second.x, DCA_sign);
0887       DCA_distance_outer_Y->Fill(i.second.y, DCA_sign);
0888     }
0889 
0890   }  // note : end of the loop for the cluster pair
0891 
0892   // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
0893   delete DCA_distance_inner_phi_peak;
0894   DCA_distance_inner_phi_peak = (TH2F*) DCA_distance_inner_phi->Clone("DCA_distance_inner_phi_peak");
0895   TH2F_threshold(DCA_distance_inner_phi_peak, 0.5);  // todo : the background cut can be modified, the ratio 0.5
0896   delete DCA_distance_inner_phi_peak_profile;
0897   DCA_distance_inner_phi_peak_profile = DCA_distance_inner_phi_peak->ProfileX("DCA_distance_inner_phi_peak_profile");
0898 
0899   double point_index = 0;
0900   std::vector<double> hist_column_content = SumTH2FColumnContent(DCA_distance_inner_phi_peak);
0901   for (int i = 0; i < DCA_distance_inner_phi_peak_profile->GetNbinsX(); i++)
0902   {
0903     if (hist_column_content[i] < 5)
0904     {
0905       continue;
0906     }  // note : in order to remove some remaining background
0907 
0908     DCA_distance_inner_phi_peak_profile_graph->SetPoint(point_index,
0909                                                         DCA_distance_inner_phi_peak_profile->GetBinCenter(i + 1),
0910                                                         DCA_distance_inner_phi_peak_profile->GetBinContent(i + 1));
0911     // std::cout<<"("<<DCA_distance_inner_phi_peak_profile->GetBinCenter(i+1)<<", "<< DCA_distance_inner_phi_peak_profile->GetBinContent(i+1)<<")"<<std::endl;
0912     point_index += 1;
0913   }
0914 
0915   //------------------------------------------------------------------
0916   // this is used to constrain the quadrant
0917   // info_vec[3];
0918   // todo : the fit range of the gaussian fit can be modified here
0919   //
0920   horizontal_fit_inner->SetParameter(0, 0);
0921   DCA_distance_inner_phi_peak_profile_graph->Fit(horizontal_fit_inner, "NQ", "", 0, 360);
0922   //------------------------------------------------------------------
0923 
0924   if (m_enable_qa)
0925   {
0926     DCA_distance_inner_phi_peak_profile_graph->Fit(gaus_fit, "NQ", "",
0927                                                    cos_fit->GetParameter(2) - guas_fit_range,   // note : what we want and need is the peak position,
0928                                                    cos_fit->GetParameter(2) + guas_fit_range);  // so we fit the peak again
0929     DCA_distance_inner_phi_peak_profile_graph->Fit(cos_fit, "NQ", "", cos_fit_rangel, cos_fit_ranger);
0930   }
0931 
0932   // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
0933   delete angle_diff_inner_phi_peak;
0934   angle_diff_inner_phi_peak = (TH2F*) angle_diff_inner_phi->Clone("angle_diff_inner_phi_peak");
0935   TH2F_threshold_advanced_2(angle_diff_inner_phi_peak, 0.5);  // todo : threshold ratio can be modified here
0936   hist_column_content = SumTH2FColumnContent(angle_diff_inner_phi_peak);
0937   angle_diff_inner_phi_peak_profile = angle_diff_inner_phi_peak->ProfileX("angle_diff_inner_phi_peak_profile");
0938   point_index = 0;
0939   for (int i = 0; i < angle_diff_inner_phi_peak_profile->GetNbinsX(); i++)
0940   {
0941     if (hist_column_content[i] < 5)
0942     {
0943       continue;
0944     }  // note : in order to remove some remaining background
0945 
0946     angle_diff_inner_phi_peak_profile_graph->SetPoint(point_index,
0947                                                       angle_diff_inner_phi_peak_profile->GetBinCenter(i + 1),
0948                                                       angle_diff_inner_phi_peak_profile->GetBinContent(i + 1));
0949     // std::cout<<"("<<angle_diff_inner_phi_peak_profile->GetBinCenter(i+1)<<", "<< angle_diff_inner_phi_peak_profile->GetBinContent(i+1)<<")"<<std::endl;
0950     point_index += 1;
0951   }
0952 
0953   //------------------------------------------------------------------
0954   // this is used to constrain the quadrant
0955   // info_vec[5];
0956   horizontal_fit_angle_diff_inner->SetParameter(0, 0);
0957   angle_diff_inner_phi_peak_profile_graph->Fit(horizontal_fit_angle_diff_inner, "NQ", "", 0, 360);
0958   //------------------------------------------------------------------
0959 
0960   angle_diff_inner_phi_peak_profile_graph->Set(0);
0961   DCA_distance_outer_phi_peak_profile_graph->Set(0);
0962 
0963   //------------------------------------------------------------------------
0964   // all others are not used for XY vertex calculation. for QA
0965   //------------------------------------------------------------------------
0966 
0967   // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
0968   delete angle_diff_new_bkg_remove;
0969   angle_diff_new_bkg_remove = (TH1F*) angle_diff_new->Clone("angle_diff_new_bkg_remove");
0970   angle_diff_new_bkg_remove->SetLineColor(2);
0971   Hist_1D_bkg_remove(angle_diff_new_bkg_remove, 1.5);
0972 
0973   if (m_enable_qa)
0974   {
0975     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
0976     delete DCA_distance_outer_phi_peak;
0977     DCA_distance_outer_phi_peak = (TH2F*) DCA_distance_outer_phi->Clone("DCA_distance_outer_phi_peak");
0978     TH2F_threshold(DCA_distance_outer_phi_peak, 0.5);  // todo : the background cut can be modified, the ratio 0.5
0979     DCA_distance_outer_phi_peak_profile = DCA_distance_outer_phi_peak->ProfileX("DCA_distance_outer_phi_peak_profile");
0980     point_index = 0;
0981     hist_column_content = SumTH2FColumnContent(DCA_distance_outer_phi_peak);
0982     for (int i = 0; i < DCA_distance_outer_phi_peak_profile->GetNbinsX(); i++)
0983     {
0984       if (hist_column_content[i] < 5)
0985       {
0986         continue;
0987       }  // note : in order to remove some remaining background
0988 
0989       DCA_distance_outer_phi_peak_profile_graph->SetPoint(point_index,
0990                                                           DCA_distance_outer_phi_peak_profile->GetBinCenter(i + 1),
0991                                                           DCA_distance_outer_phi_peak_profile->GetBinContent(i + 1));
0992       // std::cout<<"("<<DCA_distance_outer_phi_peak_profile->GetBinCenter(i+1)<<", "<< DCA_distance_outer_phi_peak_profile->GetBinContent(i+1)<<")"<<std::endl;
0993       point_index += 1;
0994     }
0995 
0996     horizontal_fit_outer->SetParameter(0, 0);
0997     // todo : the fit range of the gaussian fit can be modified here
0998     DCA_distance_outer_phi_peak_profile_graph->Fit(horizontal_fit_outer, "NQ", "", 0, 360);
0999 
1000     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1001     delete angle_diff_outer_phi_peak;
1002     angle_diff_outer_phi_peak = (TH2F*) angle_diff_outer_phi->Clone("angle_diff_outer_phi_peak");
1003     TH2F_threshold_advanced_2(angle_diff_outer_phi_peak, 0.5);  // todo : threshold ratio can be modified here
1004     hist_column_content = SumTH2FColumnContent(angle_diff_outer_phi_peak);
1005     angle_diff_outer_phi_peak_profile = angle_diff_outer_phi_peak->ProfileX("angle_diff_outer_phi_peak_profile");
1006     point_index = 0;
1007     for (int i = 0; i < angle_diff_outer_phi_peak_profile->GetNbinsX(); i++)
1008     {
1009       if (hist_column_content[i] < 5)
1010       {
1011         continue;
1012       }  // note : in order to remove some remaining background
1013 
1014       angle_diff_outer_phi_peak_profile_graph->SetPoint(point_index,
1015                                                         angle_diff_outer_phi_peak_profile->GetBinCenter(i + 1),
1016                                                         angle_diff_outer_phi_peak_profile->GetBinContent(i + 1));
1017       // std::cout<<"("<<angle_diff_outer_phi_peak_profile->GetBinCenter(i+1)<<", "<< angle_diff_outer_phi_peak_profile->GetBinContent(i+1)<<")"<<std::endl;
1018       point_index += 1;
1019     }
1020 
1021     horizontal_fit_angle_diff_outer->SetParameter(0, 0);
1022     angle_diff_outer_phi_peak_profile_graph->Fit(horizontal_fit_angle_diff_outer, "NQ", "", 0, 360);
1023 
1024     // note : ----------- ----------- ----------- ---------- ----------- ----------- ---------- ----------- ----------- ----------- -----------
1025 
1026     angle_diff_outer_phi_peak_profile_graph->Set(0);
1027     DCA_distance_inner_phi_peak_profile_graph->Set(0);
1028   }
1029 
1030   if (m_enable_qa && print_message_opt == true)
1031   {
1032     std::cout << "circle radius : " << std::abs(gaus_fit->GetParameter(0) + gaus_fit->GetParameter(3))
1033               << " possible correction angle : " << gaus_fit->GetParameter(1) << std::endl;
1034   }
1035 }
1036 
1037 // void INTTXYvtx::PrintPlotsVTXxy(std::string sub_out_folder_name)
1038 void INTTXYvtx::PrintPlotsVTXxy()
1039 {
1040   if (m_enable_drawhist)
1041   {
1042     std::string s_inttlabel = (boost::format("#it{#bf{sPHENIX INTT}} %s") % plot_text).str().c_str();
1043     std::string s_pdfname = out_folder_directory + "/" + m_quad_pdfname;
1044     std::cout << s_pdfname << std::endl;
1045 
1046     // note : -----------------------------------------------------------------------------------------
1047     DCA_distance_inner_phi->Draw("colz0");
1048     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1049     c1->Print(s_pdfname.c_str());
1050     c1->Clear();
1051 
1052     // note : -----------------------------------------------------------------------------------------
1053     angle_diff_inner_phi->Draw("colz0");
1054     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1055     c1->Print(s_pdfname.c_str());
1056     c1->Clear();
1057 
1058     // note : -----------------------------------------------------------------------------------------
1059     DCA_distance_inner_phi_peak->SetStats(false);
1060     DCA_distance_inner_phi_peak->GetXaxis()->SetTitle("Inner phi [degree]");
1061     DCA_distance_inner_phi_peak->GetYaxis()->SetTitle("DCA [mm]");
1062     DCA_distance_inner_phi_peak->Draw("colz0");
1063     DCA_distance_inner_phi_peak_profile->Draw("same");
1064     // cos_fit -> Draw("l same");
1065     // gaus_fit -> Draw("l same");
1066     horizontal_fit_inner->Draw("l same");
1067     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1068     draw_text->DrawLatex(0.25, 0.84, (boost::format("#color[2]{Assumed vertex: %.3f mm, %.3f mm}") % current_vtxX % current_vtxY).str().c_str());
1069     draw_text->DrawLatex(0.25, 0.80, (boost::format("#color[2]{Pol0 fit chi2/NDF: %.3f, fit error: %.3f}") % (horizontal_fit_inner->GetChisquare() / double(horizontal_fit_inner->GetNDF())) % horizontal_fit_inner->GetParError(0)).str().c_str());
1070     c1->Print(s_pdfname.c_str());
1071     c1->Clear();
1072 
1073     // note : -----------------------------------------------------------------------------------------
1074     angle_diff_inner_phi_peak->SetStats(false);
1075     angle_diff_inner_phi_peak->GetXaxis()->SetTitle("Inner phi [degree]");
1076     angle_diff_inner_phi_peak->GetYaxis()->SetTitle("Inner - Outer [degree]");
1077     angle_diff_inner_phi_peak->Draw("colz0");
1078     angle_diff_inner_phi_peak_profile->Draw("same");
1079     horizontal_fit_angle_diff_inner->Draw("l same");
1080     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1081     draw_text->DrawLatex(0.25, 0.84, (boost::format("#color[2]{Assumed vertex: %.3f mm, %.3f mm}") % current_vtxX % current_vtxY).str().c_str());
1082     c1->Print(s_pdfname.c_str());
1083     c1->Clear();
1084 
1085     //----------------------
1086     if (m_enable_qa)
1087     {
1088       // note : -----------------------------------------------------------------------------------------
1089       DCA_distance_outer_phi->Draw("colz0");
1090       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1091       c1->Print(s_pdfname.c_str());
1092       c1->Clear();
1093 
1094       // note : -----------------------------------------------------------------------------------------
1095       angle_diff_outer_phi->Draw("colz0");
1096       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1097       c1->Print(s_pdfname.c_str());
1098       c1->Clear();
1099 
1100       // note : -----------------------------------------------------------------------------------------
1101       angle_diff->SetMinimum(0);
1102       angle_diff->Draw("hist");
1103       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1104       c1->Print(s_pdfname.c_str());
1105       c1->Clear();
1106 
1107       // note : -----------------------------------------------------------------------------------------
1108       angle_diff_new->SetMinimum(0);
1109       angle_diff_new->Draw("hist");
1110       angle_diff_new_bkg_remove->Draw("hist same");
1111       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1112       draw_text->DrawLatex(0.4, 0.80, (boost::format("#color[2]{Dist. StdDev: %.4f}") % angle_diff_new_bkg_remove->GetStdDev()).str().c_str());
1113       c1->Print(s_pdfname.c_str());
1114       c1->Clear();
1115 
1116       // note : -----------------------------------------------------------------------------------------
1117       DCA_distance->Draw("hist");
1118       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1119       c1->Print(s_pdfname.c_str());
1120       c1->Clear();
1121 
1122       // note : -----------------------------------------------------------------------------------------
1123       angle_correlation->Draw("colz0");
1124       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1125       c1->Print(s_pdfname.c_str());
1126       c1->Clear();
1127 
1128       // note : -----------------------------------------------------------------------------------------
1129       angle_diff_DCA_dist->Draw("colz0");
1130       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1131       c1->Print(s_pdfname.c_str());
1132       c1->Clear();
1133 
1134       // note : -----------------------------------------------------------------------------------------
1135       DCA_point->Draw("colz0");
1136       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1137       c1->Print(s_pdfname.c_str());
1138       c1->Clear();
1139 
1140       // note : -----------------------------------------------------------------------------------------
1141       DCA_distance_inner_X->Draw("colz0");
1142       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1143       c1->Print(s_pdfname.c_str());
1144       c1->Clear();
1145 
1146       // note : -----------------------------------------------------------------------------------------
1147       DCA_distance_inner_Y->Draw("colz0");
1148       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1149       c1->Print(s_pdfname.c_str());
1150       c1->Clear();
1151 
1152       // note : -----------------------------------------------------------------------------------------
1153       DCA_distance_outer_X->Draw("colz0");
1154       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1155       c1->Print(s_pdfname.c_str());
1156       c1->Clear();
1157 
1158       // note : -----------------------------------------------------------------------------------------
1159       DCA_distance_outer_Y->Draw("colz0");
1160       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1161       c1->Print(s_pdfname.c_str());
1162       c1->Clear();
1163 
1164       // note : -----------------------------------------------------------------------------------------
1165       DCA_distance_outer_phi_peak->SetStats(false);
1166       DCA_distance_outer_phi_peak->GetXaxis()->SetTitle("Outer phi [degree]");
1167       DCA_distance_outer_phi_peak->GetYaxis()->SetTitle("DCA [mm]");
1168       DCA_distance_outer_phi_peak->Draw("colz0");
1169       DCA_distance_outer_phi_peak_profile->Draw("same");
1170       horizontal_fit_outer->Draw("l same");
1171       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1172       draw_text->DrawLatex(0.25, 0.80, (boost::format("#color[2]{Assumed vertex: %.3f mm, %.3f mm}") % current_vtxX % current_vtxY).str().c_str());
1173       c1->Print(s_pdfname.c_str());
1174       c1->Clear();
1175 
1176       // note : -----------------------------------------------------------------------------------------
1177       angle_diff_outer_phi_peak->SetStats(false);
1178       angle_diff_outer_phi_peak->GetXaxis()->SetTitle("Outer phi [degree]");
1179       angle_diff_outer_phi_peak->GetYaxis()->SetTitle("Inner - Outer [degree]");
1180       angle_diff_outer_phi_peak->Draw("colz0");
1181       angle_diff_outer_phi_peak_profile->Draw("same");
1182       horizontal_fit_angle_diff_outer->Draw("l same");
1183       ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, (boost::format("%s, peak : %f") % s_inttlabel.c_str() % peek).str().c_str());
1184       draw_text->DrawLatex(0.25, 0.84, (boost::format("#color[2]{Assumed vertex: %.3f mm, %.3f mm}") % current_vtxX % current_vtxY).str().c_str());
1185       c1->Print(s_pdfname.c_str());
1186       c1->Clear();
1187     }
1188   }
1189 }
1190 
1191 void INTTXYvtx::ClearHist(int /*print_option*/)
1192 {
1193   // clear histograms for quadrant method
1194   DCA_distance_inner_phi->Reset("ICESM");
1195   angle_diff_inner_phi->Reset("ICESM");
1196 
1197   if (m_enable_qa)
1198   {
1199     DCA_distance_outer_phi->Reset("ICESM");
1200     angle_diff_outer_phi->Reset("ICESM");
1201 
1202     angle_correlation->Reset("ICESM");
1203     angle_diff_DCA_dist->Reset("ICESM");
1204     angle_diff->Reset("ICESM");
1205     DCA_point->Reset("ICESM");
1206     DCA_distance->Reset("ICESM");
1207     DCA_distance_inner_X->Reset("ICESM");
1208     DCA_distance_inner_Y->Reset("ICESM");
1209     DCA_distance_outer_X->Reset("ICESM");
1210     DCA_distance_outer_Y->Reset("ICESM");
1211 
1212     angle_diff_new->Reset("ICESM");
1213     //
1214     // DCA_distance_inner_phi_peak -> Reset("ICESM");
1215     // DCA_distance_outer_phi_peak -> Reset("ICESM");
1216     // angle_diff_inner_phi_peak -> Reset("ICESM");
1217     // angle_diff_outer_phi_peak -> Reset("ICESM");
1218     // angle_diff_new_bkg_remove -> Reset("ICESM");
1219   }
1220 
1221   delete angle_diff_new_bkg_remove;
1222   angle_diff_new_bkg_remove = nullptr;
1223   delete angle_diff_inner_phi_peak;
1224   angle_diff_inner_phi_peak = nullptr;
1225   delete angle_diff_outer_phi_peak;
1226   angle_diff_outer_phi_peak = nullptr;
1227   delete DCA_distance_inner_phi_peak;
1228   DCA_distance_inner_phi_peak = nullptr;
1229   delete DCA_distance_outer_phi_peak;
1230   DCA_distance_outer_phi_peak = nullptr;
1231 }
1232 
1233 void INTTXYvtx::EndRun()
1234 {
1235   if (!m_initialized)
1236   {
1237     std::cout << "INTTXYvtx is not initialized, abort in EndRun" << std::endl;
1238     exit(1);
1239   }
1240 
1241   //------------------------
1242   // write histograms
1243   if (m_savehist)
1244   {
1245     if (!std::filesystem::exists(out_folder_directory))
1246     {
1247       std::filesystem::create_directory(out_folder_directory);
1248     }
1249 
1250     TFile* file_out = new TFile((boost::format("%s/run_XY_histo.root") % out_folder_directory).str().c_str(), "RECREATE");
1251 
1252     for (auto& itr : m_v_hist)
1253     {
1254       itr->Write();
1255     }
1256 
1257     if (xy_hist != nullptr)
1258     {
1259       xy_hist->Write();
1260     }
1261     if (xy_hist_bkgrm != nullptr)
1262     {
1263       xy_hist_bkgrm->Write();
1264     }
1265 
1266     file_out->Close();
1267     delete file_out;
1268   }
1269 
1270   // ProcessEvt QA histos
1271   //--    inner_pos_xy -> Reset("ICESM");
1272   //--    outer_pos_xy -> Reset("ICESM");
1273   //--    inner_outer_pos_xy -> Reset("ICESM");
1274   //--    N_cluster_correlation -> Reset("ICESM");
1275   //--    N_cluster_correlation_close -> Reset("ICESM");
1276 
1277   // file_out -> cd();
1278   // tree_out -> SetDirectory(file_out);
1279   // tree_out -> Write();
1280   // file_out -> Close();
1281 
1282   // c1 -> Delete();
1283   // ltx -> Delete();
1284   // draw_text -> Delete();
1285   // cos_fit -> Delete();
1286   // gaus_fit -> Delete();
1287 
1288   // horizontal_fit_inner -> Delete();
1289   // horizontal_fit_angle_diff_inner -> Delete();
1290   // horizontal_fit_outer -> Delete();
1291   // horizontal_fit_angle_diff_outer -> Delete();
1292 
1293   //--delete gROOT->FindObject("angle_diff");
1294   //--delete gROOT->FindObject("angle_diff_new");
1295 
1296   // angle_correlation -> Delete();
1297   // angle_diff_DCA_dist -> Delete();
1298   // angle_diff -> Delete();
1299   // angle_diff_new -> Delete();
1300   // angle_diff_new_bkg_remove -> Delete();
1301   // inner_pos_xy -> Delete();
1302   // outer_pos_xy -> Delete();
1303   // inner_outer_pos_xy -> Delete();
1304   // DCA_point -> Delete();
1305   // DCA_distance -> Delete();
1306   // N_cluster_correlation -> Delete();
1307   // N_cluster_correlation_close -> Delete();
1308 
1309   // DCA_distance_inner_X -> Delete();
1310   // DCA_distance_inner_Y -> Delete();
1311   // DCA_distance_outer_X -> Delete();
1312   // DCA_distance_outer_Y -> Delete();
1313 
1314   // DCA_distance_inner_phi -> Delete();
1315   // DCA_distance_inner_phi_peak -> Delete();
1316   // DCA_distance_inner_phi_peak_profile -> Delete();
1317   // DCA_distance_outer_phi -> Delete();
1318   // DCA_distance_outer_phi_peak -> Delete();
1319   // DCA_distance_outer_phi_peak_profile -> Delete();
1320 
1321   // angle_diff_inner_phi -> Delete();
1322   // angle_diff_inner_phi_peak -> Delete();
1323   // angle_diff_inner_phi_peak_profile -> Delete();
1324   // angle_diff_outer_phi -> Delete();
1325   // angle_diff_outer_phi_peak -> Delete();
1326   // angle_diff_outer_phi_peak_profile -> Delete();
1327 
1328   // DCA_distance_inner_phi_peak_final -> Delete();
1329   // angle_diff_inner_phi_peak_final -> Delete();
1330 
1331   return;
1332 }
1333 
1334 void INTTXYvtx::TH2F_threshold(TH2F* hist, double threshold)
1335 {
1336   double max_cut = hist->GetMaximum() * threshold;
1337 
1338   for (int xi = 0; xi < hist->GetNbinsX(); xi++)
1339   {
1340     for (int yi = 0; yi < hist->GetNbinsY(); yi++)
1341     {
1342       if (hist->GetBinContent(xi + 1, yi + 1) < max_cut)
1343       {
1344         hist->SetBinContent(xi + 1, yi + 1, 0);
1345       }
1346     }
1347   }
1348 }
1349 
1350 void INTTXYvtx::TH2F_threshold_advanced_2(TH2F* hist, double threshold)
1351 {
1352   // note : this function is to remove the background of the 2D histogram
1353   // note : but the threshold is given by average of the contents of the top "chosen_bin" bins and timing the threshold
1354   double max_cut = 0;
1355   int chosen_bin = 7;
1356 
1357   std::vector<float> all_bin_content_vec{};  //--all_bin_content_vec.clear();
1358   for (int xi = 0; xi < hist->GetNbinsX(); xi++)
1359   {
1360     for (int yi = 0; yi < hist->GetNbinsY(); yi++)
1361     {
1362       all_bin_content_vec.push_back(hist->GetBinContent(xi + 1, yi + 1));
1363     }
1364   }
1365   std::vector<unsigned long> ind(all_bin_content_vec.size(), 0);
1366   TMath::Sort(all_bin_content_vec.size(), &all_bin_content_vec[0], &ind[0]);
1367   for (int i = 0; i < chosen_bin; i++)
1368   {
1369     max_cut += all_bin_content_vec[ind[i]];
1370     /*std::cout<<"test : "<<all_bin_content_vec[ind[i]]<<std::endl;*/
1371   }
1372 
1373   max_cut = (max_cut / double(chosen_bin)) * threshold;
1374 
1375   for (int xi = 0; xi < hist->GetNbinsX(); xi++)
1376   {
1377     for (int yi = 0; yi < hist->GetNbinsY(); yi++)
1378     {
1379       if (hist->GetBinContent(xi + 1, yi + 1) < max_cut)
1380       {
1381         hist->SetBinContent(xi + 1, yi + 1, 0);
1382       }
1383     }
1384   }
1385 }
1386 
1387 std::vector<double>
1388 INTTXYvtx::calculateDistanceAndClosestPoint(
1389     double x1,
1390     double y1,
1391     double x2,
1392     double y2,
1393     double target_x,
1394     double target_y)
1395 {
1396   if (x1 != x2)
1397   {
1398     // Calculate the slope and intercept of the line passing through (x1, y1) and (x2, y2)
1399     double a = (y2 - y1) / (x2 - x1);
1400     double b = y1 - a * x1;
1401 
1402     // std::cout<<"slope : y="<<a<<"x+"<<b<<std::endl;
1403 
1404     // Calculate the closest distance from (target_x, target_y) to the line y = ax + b
1405     double closest_distance = std::abs(a * target_x - target_y + b) / std::sqrt(a * a + 1);
1406 
1407     // Calculate the coordinates of the closest point (Xc, Yc) on the line y = ax + b
1408     double Xc = (target_x + a * target_y - a * b) / (a * a + 1);
1409     double Yc = a * Xc + b;
1410 
1411     return {closest_distance, Xc, Yc};
1412   }
1413   else
1414   {
1415     double closest_distance = std::abs(x1 - target_x);
1416     double Xc = x1;
1417     double Yc = target_y;
1418 
1419     return {closest_distance, Xc, Yc};
1420   }
1421 }
1422 
1423 // note : Function to calculate the angle between two vectors in degrees using the cross product
1424 double INTTXYvtx::calculateAngleBetweenVectors(
1425     double x1,
1426     double y1,
1427     double x2,
1428     double y2,
1429     double targetX,
1430     double targetY)
1431 {
1432   // Calculate the vectors vector_1 (point_1 to point_2) and vector_2 (point_1 to target)
1433   double vector1X = x2 - x1;
1434   double vector1Y = y2 - y1;
1435 
1436   double vector2X = targetX - x1;
1437   double vector2Y = targetY - y1;
1438 
1439   // Calculate the cross product of vector_1 and vector_2 (z-component)
1440   double crossProduct = vector1X * vector2Y - vector1Y * vector2X;
1441 
1442   // std::cout<<" crossProduct : "<<crossProduct<<std::endl;
1443 
1444   // Calculate the magnitudes of vector_1 and vector_2
1445   double magnitude1 = std::sqrt(vector1X * vector1X + vector1Y * vector1Y);
1446   double magnitude2 = std::sqrt(vector2X * vector2X + vector2Y * vector2Y);
1447 
1448   // Calculate the angle in radians using the inverse tangent of the cross product and dot product
1449   //--double dotProduct = vector1X * vector2X + vector1Y * vector2Y;
1450 
1451   //--double angleInRadians = std::atan2(std::abs(crossProduct), dotProduct);
1452   // Convert the angle from radians to degrees and return it
1453   //--double angleInDegrees = angleInRadians * 180.0 / M_PI;
1454 
1455   double angleInRadians_new = std::asin(crossProduct / (magnitude1 * magnitude2));
1456 
1457   //--double angleInDegrees_new = angleInRadians_new * 180.0 / M_PI;
1458 
1459   // std::cout<<"angle : "<<angleInDegrees_new<<std::endl;
1460 
1461   double DCA_value = sin(angleInRadians_new) * magnitude2;  // DCA_value insread of DCA_distance
1462 
1463   return DCA_value;
1464 }
1465 
1466 void INTTXYvtx::Hist_1D_bkg_remove(TH1F* hist_in, double factor)
1467 {
1468   // todo : N bins considered to be used in the background quantification
1469   std::vector<double> Nbin_content_vec{};
1470   for (int i = hist_in->GetNbinsX() - 5; i < hist_in->GetNbinsX(); i++)
1471   {
1472     Nbin_content_vec.push_back(hist_in->GetBinContent(i + 1));
1473   }
1474 
1475   double bkg_level = accumulate(Nbin_content_vec.begin(), Nbin_content_vec.end(), 0.0) / Nbin_content_vec.size();
1476   // std::cout<<"test, bkg cut : "<<bkg_level * factor<<std::endl;
1477 
1478   for (int i = 0; i < hist_in->GetNbinsX(); i++)
1479   {
1480     // note : the background rejection is here : bkg_level * 1.5 for the time being
1481     double bin_content = (hist_in->GetBinContent(i + 1) <= bkg_level * factor)
1482                              ? 0.
1483                              : (hist_in->GetBinContent(i + 1) - bkg_level * factor);
1484 
1485     hist_in->SetBinContent(i + 1, bin_content);
1486   }
1487 }
1488 
1489 void INTTXYvtx::DrawTGraphErrors(
1490     std::vector<double> x_vec,
1491     std::vector<double> y_vec,
1492     std::vector<double> xE_vec,
1493     std::vector<double> yE_vec,
1494     const std::string& output_directory,
1495     std::vector<std::string> plot_name)
1496 {
1497   if (m_enable_drawhist)
1498   {
1499     c1->cd();
1500 
1501     TGraphErrors* g = new TGraphErrors(x_vec.size(), &x_vec[0], &y_vec[0], &xE_vec[0], &yE_vec[0]);
1502     g->SetMarkerStyle(20);
1503     g->SetMarkerSize(1.5);
1504     g->SetMarkerColor(1);
1505     g->SetLineColor(1);
1506     g->GetXaxis()->SetTitle(plot_name[1].c_str());
1507     g->GetYaxis()->SetTitle(plot_name[2].c_str());
1508     g->SetTitle(plot_name[0].c_str());
1509     if (plot_name.size() == 4)
1510     {
1511       g->Draw(plot_name[3].c_str());
1512     }
1513     else
1514     {
1515       g->Draw("AP");
1516     }
1517 
1518     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, (boost::format("#it{#bf{sPHENIX INTT}} %s") % plot_text).str().c_str());
1519     c1->Print((boost::format("%s/%s") % output_directory.c_str() % plot_name[4]).str().c_str());
1520     c1->Clear();
1521 
1522     delete g;
1523   }
1524 }
1525 
1526 void INTTXYvtx::Draw2TGraph(
1527     std::vector<double> x1_vec,
1528     std::vector<double> y1_vec,
1529     std::vector<double> x2_vec,
1530     std::vector<double> y2_vec,
1531     const std::string& output_directory,
1532     std::vector<std::string> plot_name)
1533 {
1534   if (m_enable_drawhist)
1535   {
1536     c1->cd();
1537     c1->SetLogy(1);
1538 
1539     TGraph* g1 = new TGraph(x1_vec.size(), &x1_vec[0], &y1_vec[0]);
1540     g1->SetMarkerStyle(5);
1541     g1->SetMarkerSize(1);
1542     g1->SetMarkerColor(1);
1543     g1->SetLineColor(1);
1544     g1->GetXaxis()->SetTitle(plot_name[1].c_str());
1545     g1->GetYaxis()->SetTitle(plot_name[2].c_str());
1546     g1->GetXaxis()->SetNdivisions(505);
1547     g1->GetXaxis()->SetLimits(-1, x1_vec[x1_vec.size() - 1] + 1);
1548     g1->SetTitle(plot_name[0].c_str());
1549     g1->Draw("AP");
1550 
1551     TGraph* g2 = new TGraph(x2_vec.size(), &x2_vec[0], &y2_vec[0]);
1552     g2->SetMarkerStyle(5);
1553     g2->SetMarkerSize(1);
1554     g2->SetMarkerColor(2);
1555     g2->SetLineColor(2);
1556     g2->Draw("PL same");
1557 
1558     TLegend* legend = new TLegend(0.4, 0.75, 0.65, 0.9);
1559     // legend -> SetMargin(0);
1560     legend->SetTextSize(0.03);
1561     legend->AddEntry(g1, "Tested vertex candidates", "p");
1562     legend->AddEntry(g2, "Better performed candidates", "p");
1563     legend->Draw("same");
1564     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, (boost::format("#it{#bf{sPHENIX INTT}} %s") % plot_text).str().c_str());
1565     c1->Print((boost::format("%s/%s") % output_directory.c_str() % plot_name[4]).str().c_str());
1566     c1->Clear();
1567     c1->SetLogy(0);
1568 
1569     delete g1;
1570     delete g2;
1571     delete legend;
1572   }
1573 }
1574 
1575 std::vector<double> INTTXYvtx::SumTH2FColumnContent(TH2F* hist_in)
1576 {
1577   std::vector<double> sum_vec;
1578   sum_vec.clear();
1579   for (int i = 0; i < hist_in->GetNbinsX(); i++)
1580   {
1581     double sum = 0;
1582     for (int j = 0; j < hist_in->GetNbinsY(); j++)
1583     {
1584       sum += hist_in->GetBinContent(i + 1, j + 1);
1585     }
1586     sum_vec.push_back(sum);
1587   }
1588   return sum_vec;
1589 }
1590 
1591 std::vector<std::pair<double, double>> INTTXYvtx::Get4vtx(std::pair<double, double> origin, double length)
1592 {
1593   std::vector<std::pair<double, double>> unit_vtx = {{1, 1}, {-1, 1}, {-1, -1}, {1, -1}};
1594   std::vector<std::pair<double, double>> vec_out{};  //-- vec_out.clear();
1595 
1596   vec_out.reserve(unit_vtx.size());
1597   for (std::pair i1 : unit_vtx)
1598   {
1599     vec_out.emplace_back(i1.first * length + origin.first, i1.second * length + origin.second);
1600   }
1601 
1602   return vec_out;
1603 }
1604 
1605 void INTTXYvtx::TH2F_FakeClone(TH2F* hist_in, TH2F* hist_out)
1606 {
1607   if (hist_in->GetNbinsX() != hist_out->GetNbinsX() ||
1608       hist_in->GetNbinsY() != hist_out->GetNbinsY())
1609   {
1610     std::cout << "In INTTXYvtx::TH2F_FakeClone, the input and output histogram have different binning!" << std::endl;
1611     return;
1612   }
1613 
1614   for (int i = 0; i < hist_in->GetNbinsX(); i++)
1615   {
1616     for (int j = 0; j < hist_in->GetNbinsY(); j++)
1617     {
1618       hist_out->SetBinContent(i + 1, j + 1, hist_in->GetBinContent(i + 1, j + 1));
1619     }
1620   }
1621 }
1622 
1623 void INTTXYvtx::TH1F_FakeClone(TH1F* hist_in, TH1F* hist_out)
1624 {
1625   if (hist_in->GetNbinsX() != hist_out->GetNbinsX())
1626   {
1627     std::cout << "In INTTXYvtx::TH1F_FakeClone, the input and output histogram have different binning!" << std::endl;
1628     return;
1629   }
1630 
1631   for (int i = 0; i < hist_in->GetNbinsX(); i++)
1632   {
1633     hist_out->SetBinContent(i + 1, hist_in->GetBinContent(i + 1));
1634   }
1635 }
1636 
1637 void INTTXYvtx::TH2FSampleLineFill(
1638     TH2F* hist_in,
1639     double segmentation,
1640     std::pair<double, double> inner_clu,
1641     std::pair<double, double> outer_clu)
1642 {
1643   if (!m_initialized)
1644   {
1645     std::cout << "INTTXYvtx is not initialized, abort in MacroVTXSquare" << std::endl;
1646     exit(1);
1647   }
1648 
1649   double x_min = hist_in->GetXaxis()->GetXmin();
1650   double x_max = hist_in->GetXaxis()->GetXmax();
1651   double y_min = hist_in->GetYaxis()->GetXmin();
1652   double y_max = hist_in->GetYaxis()->GetXmax();
1653 
1654   double seg_x, seg_y;
1655   double angle;
1656   int n_seg = 0;
1657 
1658   while (true)
1659   {
1660     angle = atan2(inner_clu.second - outer_clu.second, inner_clu.first - outer_clu.first);
1661     seg_x = (n_seg * segmentation) * cos(angle) + outer_clu.first;  // note : atan2(y,x), point.first is the radius
1662     seg_y = (n_seg * segmentation) * sin(angle) + outer_clu.second;
1663 
1664     if ((seg_x > x_min && seg_x < x_max && seg_y > y_min && seg_y < y_max) != true)
1665     {
1666       break;
1667     }
1668 
1669     hist_in->Fill(seg_x, seg_y);
1670     n_seg += 1;
1671   }
1672 
1673   n_seg = 1;
1674   while (true)
1675   {
1676     angle = atan2(inner_clu.second - outer_clu.second, inner_clu.first - outer_clu.first);
1677     seg_x = (-1 * n_seg * segmentation) * cos(angle) + outer_clu.first;  // note : atan2(y,x), point.first is the radius
1678     seg_y = (-1 * n_seg * segmentation) * sin(angle) + outer_clu.second;
1679 
1680     if ((seg_x > x_min && seg_x < x_max && seg_y > y_min && seg_y < y_max) != true)
1681     {
1682       break;
1683     }
1684     hist_in->Fill(seg_x, seg_y);
1685     n_seg += 1;
1686   }
1687 }
1688 
1689 std::vector<std::pair<double, double>>
1690 INTTXYvtx::FillLine_FindVertex(
1691     std::pair<double, double> window_center,
1692     double segmentation,
1693     double window_width,
1694     int N_bins)
1695 {
1696   bool draw_plot = m_enable_drawhist;
1697 
1698   if (cluster_pair_vec.size() == 0)
1699   {  // minimum tracklet cut. should be tuned
1700     return {beam_origin,
1701             {0, 0},
1702             {0, 0}};
1703   }
1704 
1705   delete xy_hist;  // if xy_hist is nullptr, nothing happen;
1706   xy_hist = new TH2F("xy_hist", "xy_hist",
1707                      N_bins, -1 * window_width / 2. + window_center.first,
1708                      window_width / 2. + window_center.first,
1709                      N_bins, -1 * window_width / 2. + window_center.second,
1710                      window_width / 2. + window_center.second);
1711   xy_hist->SetStats(false);
1712   xy_hist->GetXaxis()->SetTitle("X axis [mm]");
1713   xy_hist->GetYaxis()->SetTitle("Y axis [mm]");
1714   xy_hist->GetXaxis()->SetNdivisions(505);
1715 
1716   delete xy_hist_bkgrm;
1717   xy_hist_bkgrm = new TH2F("xy_hist_bkgrm", "xy_hist_bkgrm",
1718                            N_bins, -1 * window_width / 2. + window_center.first,
1719                            window_width / 2. + window_center.first,
1720                            N_bins, -1 * window_width / 2. + window_center.second,
1721                            window_width / 2. + window_center.second);
1722   xy_hist_bkgrm->SetStats(false);
1723   xy_hist_bkgrm->GetXaxis()->SetTitle("X axis [mm]");
1724   xy_hist_bkgrm->GetYaxis()->SetTitle("Y axis [mm]");
1725   xy_hist_bkgrm->GetXaxis()->SetNdivisions(505);
1726 
1727   // std::cout<<"test test size and bin of the hist xy_hist : "<<xy_hist -> GetNbinsX()<<" "<<xy_hist -> GetNbinsY()<<std::endl;
1728   // std::cout<<"test test bin width of the hist xy_hist : "<<xy_hist -> GetXaxis() -> GetBinWidth(1)<<" "<<xy_hist -> GetYaxis() -> GetBinWidth(1)<<std::endl;
1729   // std::cout<<"draw_plot status : "<<draw_plot<<std::endl;
1730 
1731   for (auto& i : cluster_pair_vec)
1732   {
1733     std::vector<double> DCA_info_vec = calculateDistanceAndClosestPoint(
1734         i.first.x, i.first.y,
1735         i.second.x, i.second.y,
1736         window_center.first, window_center.second);
1737 
1738     double DCA_sign = calculateAngleBetweenVectors(
1739         i.second.x, i.second.y,
1740         i.first.x, i.first.y,
1741         window_center.first, window_center.second);
1742 
1743     if (DCA_info_vec[0] != fabs(DCA_sign) && fabs(DCA_info_vec[0] - fabs(DCA_sign)) > 0.1)
1744     {
1745       std::cout << "different DCA : " << DCA_info_vec[0] << " " << DCA_sign << " diff : " << DCA_info_vec[0] - fabs(DCA_sign) << std::endl;
1746     }
1747 
1748     Clus_InnerPhi_Offset = (i.first.y - window_center.second < 0)
1749                                ? atan2(i.first.y - window_center.second,
1750                                        i.first.x - window_center.first) *
1751                                          (180. / M_PI) +
1752                                      360
1753                                : atan2(i.first.y - window_center.second,
1754                                        i.first.x - window_center.first) *
1755                                      (180. / M_PI);
1756     Clus_OuterPhi_Offset = (i.second.y - window_center.second < 0)
1757                                ? atan2(i.second.y - window_center.second,
1758                                        i.second.x - window_center.first) *
1759                                          (180. / M_PI) +
1760                                      360
1761                                : atan2(i.second.y - window_center.second,
1762                                        i.second.x - window_center.first) *
1763                                      (180. / M_PI);
1764 
1765     if (fabs(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset) < 5)
1766     {
1767       TH2FSampleLineFill(xy_hist,
1768                          segmentation,
1769                          {i.first.x, i.first.y},
1770                          {DCA_info_vec[1], DCA_info_vec[2]});
1771       // note : the DCA cut may be biased since the DCA was calculated with respect to the vertex calculated by the quadrant method
1772       // if (DCA_cut.first < DCA_sign && DCA_sign < DCA_cut.second)
1773       // {
1774       //     TH2FSampleLineFill(xy_hist,
1775       //                        segmentation,
1776       //                        {cluster_pair_vec[i].first.x, cluster_pair_vec[i].first.y},
1777       //                        {DCA_info_vec[1], DCA_info_vec[2]}
1778       //                       );
1779       // }
1780     }
1781   }
1782 
1783   // note : try to remove the background
1784   TH2F_FakeClone(xy_hist, xy_hist_bkgrm);
1785   TH2F_threshold_advanced_2(xy_hist_bkgrm, 0.7);
1786 
1787   double reco_vtx_x = xy_hist_bkgrm->GetMean(1);  // note : the TH2F calculate the GetMean based on the bin center, no need to apply additional offset
1788   double reco_vtx_y = xy_hist_bkgrm->GetMean(2);  // note : the TH2F calculate the GetMean based on the bin center, no need to apply additional offset
1789 
1790   // std::cout<<"test : in the line filled, the process is almost done"<<std::endl;
1791 
1792   if (draw_plot)
1793   {
1794     TGraph* reco_vertex_gr = new TGraph();
1795     reco_vertex_gr->SetMarkerStyle(50);
1796     reco_vertex_gr->SetMarkerColor(2);
1797     reco_vertex_gr->SetMarkerSize(1);
1798     reco_vertex_gr->SetPoint(reco_vertex_gr->GetN(), reco_vtx_x, reco_vtx_y);
1799 
1800     std::string s_inttlabel = (boost::format("#it{#bf{sPHENIX INTT}} %s") % plot_text).str().c_str();
1801     // note : -----------------------------------------------------------------------------------------
1802     xy_hist->Draw("colz0");
1803     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1804     c1->Print((boost::format("%s/linefill_qa.pdf(") % out_folder_directory).str().c_str());
1805     c1->Clear();
1806 
1807     // note : -----------------------------------------------------------------------------------------
1808     xy_hist_bkgrm->Draw("colz0");
1809     draw_text->DrawLatex(0.21, 0.71 + 0.13, (boost::format("Vertex of the Run: %.4f mm, %.4f mm") % reco_vtx_x % reco_vtx_y).str().c_str());
1810     draw_text->DrawLatex(0.21, 0.67 + 0.13, (boost::format("Vertex error: %.4f mm, %.4f mm") % xy_hist_bkgrm->GetMeanError(1) % xy_hist_bkgrm->GetMeanError(2)).str().c_str());
1811     reco_vertex_gr->Draw("p same");
1812     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, s_inttlabel.c_str());
1813     c1->Print((boost::format("%s/linefill_qa.pdf)") % out_folder_directory).str().c_str());
1814     c1->Clear();
1815 
1816     // std::cout<<"test : hello, can you see me ?"<<std::endl;
1817   }
1818 
1819   return {{reco_vtx_x, reco_vtx_y},
1820           {xy_hist_bkgrm->GetMeanError(1), xy_hist_bkgrm->GetMeanError(2)},
1821           {xy_hist_bkgrm->GetStdDev(1), xy_hist_bkgrm->GetStdDev(2)}};
1822 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!