Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:12:32 

0001 #include "AvgVtxXY.h"
0002 
0003 AvgVtxXY::AvgVtxXY(
0004     int process_id_in,
0005     int runnumber_in,
0006     int run_nEvents_in,
0007     std::string input_directory_in,
0008     std::string input_file_name_in,
0009     std::string output_directory_in,
0010     std::string output_file_name_suffix_in,
0011 
0012     // note : for the event selection
0013     std::pair<double,double> MBD_vtxZ_cut_in,
0014     std::pair<int,int> INTTNClus_cut_in,
0015     int ClusAdc_cut_in,
0016     int ClusPhiSize_cut_in,
0017 
0018     bool HaveGeoOffsetTag_in,
0019     double random_range_XYZ_in,
0020     int random_seed_in,
0021     std::string input_offset_map_in // note : full map
0022 ):
0023     process_id(process_id_in),
0024     runnumber(runnumber_in),
0025     run_nEvents(run_nEvents_in),
0026     input_directory(input_directory_in),
0027     input_file_name(input_file_name_in),
0028     output_directory(output_directory_in),
0029     output_file_name_suffix(output_file_name_suffix_in),
0030     MBD_vtxZ_cut(MBD_vtxZ_cut_in),
0031     INTTNClus_cut(INTTNClus_cut_in),
0032     ClusAdc_cut(ClusAdc_cut_in),
0033     ClusPhiSize_cut(ClusPhiSize_cut_in),
0034     HaveGeoOffsetTag(HaveGeoOffsetTag_in),
0035     random_range_XYZ(random_range_XYZ_in),
0036     random_seed(random_seed_in),
0037     input_offset_map(input_offset_map_in),
0038 
0039     hopeless_pair_angle_cut(45),                           // todo: // note : degree
0040     rough_DeltaPhi_cut_degree(7),                          // todo: 
0041     InnerOuterPhi_DeltaPhi_DCA_Threshold(0.5),             // todo:
0042     LineFill_Hist_Threshold(0.7),                          // todo:
0043     InnerPhi_DeltaPhi_DCA_Fit_Tgr_HitContent_Threshold(5), // todo:
0044     TH2D_threshold_advanced_2_CheckNBins(7)                // todo:
0045 {
0046     ReadRootFile();
0047     run_nEvents = (run_nEvents == -1) ? tree_in->GetEntries() : run_nEvents;
0048     run_nEvents = (run_nEvents > tree_in->GetEntries()) ? tree_in->GetEntries() : run_nEvents;
0049     
0050     InitOutRootFile();   
0051 
0052     fit_innerphi_DeltaPhi = new TF1("fit_innerphi_DeltaPhi","pol0",-M_PI, M_PI);
0053     fit_innerphi_DCA = new TF1("fit_innerphi_DCA","pol0",-M_PI, M_PI);
0054 
0055     h2D_map.clear();
0056     h1D_map.clear();
0057     gr_map.clear();
0058     profile_map.clear();
0059 
0060     for (int layer_i = B0L0_index; layer_i <= B1L1_index; layer_i++){
0061         double N_layer_ladder = (layer_i == B0L0_index || layer_i == B0L1_index) ? nLadder_inner : nLadder_outer;
0062         for (int phi_i = 0; phi_i < N_layer_ladder; phi_i++){
0063             geo_offset_map[Form("%i_%i", layer_i, phi_i)] = {0, 0, 0};
0064         }
0065     }
0066 
0067     if (HaveGeoOffsetTag) {
0068         if (input_offset_map == "no_map"){
0069             std::cout<<"HaveGeoOffsetTag: "<<HaveGeoOffsetTag<<", but no_map from the input, therefore, generate one"<<std::endl;
0070             GenGeoOffset();
0071         }
0072         else {
0073             std::cout<<"HaveGeoOffsetTag: "<<HaveGeoOffsetTag<<", with the map: "<<input_offset_map<<std::endl;
0074             SetGeoOffset();
0075         }
0076     }
0077 
0078     // todo: it works if all the events in the files are used
0079     // out_start_evt = process_id * run_nEvents;
0080     // out_end_evt   = out_start_evt + run_nEvents - 1;
0081     
0082 }
0083 
0084 void AvgVtxXY::ReadRootFile()
0085 {
0086     file_in = TFile::Open((input_directory + "/" + input_file_name).c_str());
0087 
0088     if (file_in == nullptr) {
0089         std::cout << "Error : cannot open the file : " << input_directory + "/" + input_file_name << std::endl;
0090         exit(1);
0091     }
0092 
0093     tree_in = (TTree*)file_in->Get(input_tree_name.c_str());
0094 
0095     tree_in -> SetBranchStatus("*", 0);
0096     tree_in -> SetBranchStatus("MBD_z_vtx", 1);
0097     tree_in -> SetBranchStatus("is_min_bias", 1);
0098     tree_in -> SetBranchStatus("MBD_centrality", 1);
0099     if (runnumber != -1) {
0100         tree_in -> SetBranchStatus("InttBco_IsToBeRemoved", 1);
0101     }
0102 
0103     tree_in -> SetBranchStatus("NClus", 1);
0104     tree_in -> SetBranchStatus("ClusX", 1);
0105     tree_in -> SetBranchStatus("ClusY", 1);
0106     tree_in -> SetBranchStatus("ClusLayer", 1);
0107     tree_in -> SetBranchStatus("ClusLadderPhiId", 1);
0108     tree_in -> SetBranchStatus("ClusAdc", 1);
0109     tree_in -> SetBranchStatus("ClusPhiSize", 1);
0110 
0111     ClusX = 0;
0112     ClusY = 0;
0113     ClusLayer = 0;
0114     ClusLadderPhiId = 0;
0115     ClusAdc = 0;
0116     ClusPhiSize = 0;
0117 
0118     tree_in -> SetBranchAddress("MBD_z_vtx", &MBD_z_vtx);
0119     tree_in -> SetBranchAddress("is_min_bias", &is_min_bias);
0120     tree_in -> SetBranchAddress("MBD_centrality", &MBD_centrality);
0121     if (runnumber != -1) {
0122         tree_in -> SetBranchAddress("InttBco_IsToBeRemoved", &InttBco_IsToBeRemoved);
0123     }
0124 
0125     tree_in -> SetBranchAddress("NClus", &NClus);
0126     tree_in -> SetBranchAddress("ClusX", &ClusX);
0127     tree_in -> SetBranchAddress("ClusY", &ClusY);
0128     tree_in -> SetBranchAddress("ClusLayer", &ClusLayer);
0129     tree_in -> SetBranchAddress("ClusLadderPhiId", &ClusLadderPhiId);
0130     tree_in -> SetBranchAddress("ClusAdc", &ClusAdc);
0131     tree_in -> SetBranchAddress("ClusPhiSize", &ClusPhiSize);
0132 }
0133 
0134 void AvgVtxXY::InitOutRootFile()
0135 {
0136     std::string job_index = std::to_string( process_id );
0137     int job_index_len = 5;
0138     job_index.insert(0, job_index_len - job_index.size(), '0');
0139 
0140     std::string runnumber_str = std::to_string( runnumber );
0141     if (runnumber != -1){
0142         int runnumber_str_len = 8;
0143         runnumber_str.insert(0, runnumber_str_len - runnumber_str.size(), '0');
0144     }
0145 
0146     if (output_file_name_suffix.size() > 0 && output_file_name_suffix[0] != '_') {
0147         output_file_name_suffix = "_" + output_file_name_suffix;
0148     }
0149 
0150     output_filename = "AvgVtxXY";
0151     output_filename = (runnumber != -1) ? "Data_" + output_filename : "MC_" + output_filename;
0152     output_filename += (HaveGeoOffsetTag) ? "_GeoOffset" : "";
0153     output_filename += output_file_name_suffix;
0154     output_filename += (runnumber != -1) ? Form("_%s_%s.root",runnumber_str.c_str(),job_index.c_str()) : Form("_%s.root",job_index.c_str());
0155 
0156     file_out = new TFile((output_directory + "/" + output_filename).c_str(), "RECREATE");
0157     
0158     // note : this is for the final
0159     tree_vtxXY = new TTree("tree_vtxXY", "tree_vtxXY");
0160     tree_vtxXY -> Branch("vtxX_quadrant", &out_vtxX_quadrant);
0161     tree_vtxXY -> Branch("vtxXE_quadrant", &out_vtxXE_quadrant);
0162     tree_vtxXY -> Branch("vtxY_quadrant", &out_vtxY_quadrant);
0163     tree_vtxXY -> Branch("vtxYE_quadrant", &out_vtxYE_quadrant);
0164 
0165     tree_vtxXY -> Branch("vtxX_LineFill", &out_vtxX_LineFill);
0166     tree_vtxXY -> Branch("vtxXE_LineFill", &out_vtxXE_LineFill);
0167     tree_vtxXY -> Branch("vtxXStdDev_LineFill", &out_vtxXStdDev_LineFill);
0168     tree_vtxXY -> Branch("vtxY_LineFill", &out_vtxY_LineFill);
0169     tree_vtxXY -> Branch("vtxYE_LineFill", &out_vtxYE_LineFill);
0170     tree_vtxXY -> Branch("vtxYStdDev_LineFill", &out_vtxYStdDev_LineFill);
0171     tree_vtxXY -> Branch("run_nEvents", &out_run_nEvents);
0172     tree_vtxXY -> Branch("job_index", &out_job_index);
0173     tree_vtxXY -> Branch("file_total_event", &out_file_total_event);
0174     // tree_vtxXY -> Branch("start_evt", &out_start_evt);
0175     // tree_vtxXY -> Branch("end_evt", &out_end_evt);
0176 
0177     out_job_index = process_id; 
0178     out_file_total_event = tree_in->GetEntries();
0179     
0180     // note : for the geo offset
0181     if (HaveGeoOffsetTag) {
0182 
0183         tree_geooffset = new TTree("tree_geooffset", "tree_geooffset");
0184         tree_geooffset -> Branch("LayerID", &out_LayerID);
0185         tree_geooffset -> Branch("LadderPhiID", &out_LadderPhiID);
0186         tree_geooffset -> Branch("offset_x", &out_offset_x);
0187         tree_geooffset -> Branch("offset_y", &out_offset_y);   
0188         tree_geooffset -> Branch("offset_z", &out_offset_z);
0189     }
0190 
0191     // note : for the quadrant detail
0192     tree_quadrant_detail = new TTree("tree_quadrant_detail", "tree_quadrant_detail");
0193     tree_quadrant_detail -> Branch("iteration", &out_iteration);
0194     tree_quadrant_detail -> Branch("quadrant", &out_quadrant);
0195     tree_quadrant_detail -> Branch("assume_center_x", &out_assume_center_x);
0196     tree_quadrant_detail -> Branch("assume_center_y", &out_assume_center_y);
0197     tree_quadrant_detail -> Branch("Fit_InnerPhi_DeltaPhi_RChi2", &out_Fit_InnerPhi_DeltaPhi_RChi2);
0198     tree_quadrant_detail -> Branch("Fit_InnerPhi_DeltaPhi_Err0", &out_Fit_InnerPhi_DeltaPhi_Err0);
0199     tree_quadrant_detail -> Branch("Fit_InnerPhi_DeltaPhi_LineY", &out_Fit_InnerPhi_DeltaPhi_LineY);
0200     tree_quadrant_detail -> Branch("Fit_InnerPhi_DCA_RChi2", &out_Fit_InnerPhi_DCA_RChi2);
0201     tree_quadrant_detail -> Branch("Fit_InnerPhi_DCA_Err0", &out_Fit_InnerPhi_DCA_Err0);
0202     tree_quadrant_detail -> Branch("Fit_InnerPhi_DCA_LineY", &out_Fit_InnerPhi_DCA_LineY);
0203 
0204     // note : for the output paramteres
0205     tree_parameters = new TTree("tree_parameters", "tree_parameters");
0206     out_run_nEvents = run_nEvents;
0207     out_MBD_vtxZ_cut_l = MBD_vtxZ_cut.first;
0208     out_MBD_vtxZ_cut_r = MBD_vtxZ_cut.second;
0209     out_INTTNClus_cut_l = INTTNClus_cut.first;
0210     out_INTTNClus_cut_r = INTTNClus_cut.second;
0211     out_ClusAdc_cut = ClusAdc_cut;
0212     out_ClusPhiSize_cut = ClusPhiSize_cut;
0213     out_HaveGeoOffsetTag = HaveGeoOffsetTag;
0214     out_random_range_XYZ = random_range_XYZ;
0215     out_random_seed = random_seed;
0216     out_input_offset_map = (input_offset_map == "no_map") ? 0 : 1;
0217     out_hopeless_pair_angle_cut = hopeless_pair_angle_cut;
0218     out_rough_DeltaPhi_cut_degree = rough_DeltaPhi_cut_degree;
0219     out_InnerOuterPhi_DeltaPhi_DCA_Threshold = InnerOuterPhi_DeltaPhi_DCA_Threshold;
0220     out_LineFill_Hist_Threshold = LineFill_Hist_Threshold;
0221     out_InnerPhi_DeltaPhi_DCA_Fit_Tgr_HitContent_Threshold = InnerPhi_DeltaPhi_DCA_Fit_Tgr_HitContent_Threshold;
0222     out_TH2D_threshold_advanced_2_CheckNBins = TH2D_threshold_advanced_2_CheckNBins;
0223 
0224     tree_parameters -> Branch("run_nEvents", &out_run_nEvents);
0225     tree_parameters -> Branch("MBD_vtxZ_cut_l", &out_MBD_vtxZ_cut_l);
0226     tree_parameters -> Branch("MBD_vtxZ_cut_r", &out_MBD_vtxZ_cut_r);
0227     tree_parameters -> Branch("INTTNClus_cut_l", &out_INTTNClus_cut_l);
0228     tree_parameters -> Branch("INTTNClus_cut_r", &out_INTTNClus_cut_r);
0229     tree_parameters -> Branch("ClusAdc_cut", &out_ClusAdc_cut);
0230     tree_parameters -> Branch("ClusPhiSize_cut", &out_ClusPhiSize_cut);
0231     tree_parameters -> Branch("HaveGeoOffsetTag", &out_HaveGeoOffsetTag);
0232     tree_parameters -> Branch("random_range_XYZ", &out_random_range_XYZ);
0233     tree_parameters -> Branch("random_seed", &out_random_seed);
0234     tree_parameters -> Branch("input_offset_map", &out_input_offset_map);
0235     tree_parameters -> Branch("hopeless_pair_angle_cut", &out_hopeless_pair_angle_cut);
0236     tree_parameters -> Branch("rough_DeltaPhi_cut_degree", &out_rough_DeltaPhi_cut_degree);
0237     tree_parameters -> Branch("InnerOuterPhi_DeltaPhi_DCA_Threshold", &out_InnerOuterPhi_DeltaPhi_DCA_Threshold);
0238     tree_parameters -> Branch("LineFill_Hist_Threshold", &out_LineFill_Hist_Threshold);
0239     tree_parameters -> Branch("InnerPhi_DeltaPhi_DCA_Fit_Tgr_HitContent_Threshold", &out_InnerPhi_DeltaPhi_DCA_Fit_Tgr_HitContent_Threshold);
0240     tree_parameters -> Branch("TH2D_threshold_advanced_2_CheckNBins", &out_TH2D_threshold_advanced_2_CheckNBins);
0241 
0242     tree_parameters -> Fill();
0243 }
0244 
0245 void AvgVtxXY::geo_offset_correction_fill()
0246 {
0247     // note : it's possible that the whole system has a systematic offset, so we need to correct it
0248     // std::vector<double> XYZ_correction = offset_correction(geo_offset_map);
0249     // for (const auto& pair : geo_offset_map)
0250     // {
0251     //     geo_offset_map[pair.first] = {pair.second[0] - XYZ_correction[0], pair.second[1] - XYZ_correction[1], pair.second[2] - XYZ_correction[2]};
0252     // }
0253 
0254     for (const auto& pair : geo_offset_map)
0255     {
0256         out_LayerID = std::stoi(pair.first.substr(0,pair.first.find("_")));
0257         out_LadderPhiID = std::stoi(pair.first.substr(pair.first.find("_")+1));
0258         out_offset_x = pair.second[0];
0259         out_offset_y = pair.second[1];
0260         out_offset_z = pair.second[2];
0261         tree_geooffset -> Fill();
0262     }
0263 }
0264 
0265 void AvgVtxXY::SetGeoOffset()
0266 {
0267     geo_offset_map.clear();
0268     TFile * file_offset = TFile::Open((input_offset_map).c_str());
0269     TTree * tree_offset = (TTree*)file_offset->Get("tree_geooffset");
0270     int in_LayerID;
0271     int in_LadderPhiID;
0272     double in_offset_x;
0273     double in_offset_y;
0274     double in_offset_z;
0275 
0276     tree_offset -> SetBranchAddress("LayerID", &in_LayerID);
0277     tree_offset -> SetBranchAddress("LadderPhiID", &in_LadderPhiID);
0278     tree_offset -> SetBranchAddress("offset_x", &in_offset_x);
0279     tree_offset -> SetBranchAddress("offset_y", &in_offset_y);
0280     tree_offset -> SetBranchAddress("offset_z", &in_offset_z);
0281 
0282     for (int i = 0; i < tree_offset->GetEntries(); i++)
0283     {
0284         tree_offset -> GetEntry(i);
0285         geo_offset_map[Form("%i_%i", in_LayerID, in_LadderPhiID)] = {in_offset_x, in_offset_y, in_offset_z};
0286     }
0287 
0288     geo_offset_correction_fill();
0289 
0290     file_offset -> Close();
0291 }
0292 
0293 void AvgVtxXY::GenGeoOffset()
0294 {
0295     geo_rand = new TRandom3(0);
0296     if (random_seed != -999) {
0297         geo_rand -> SetSeed(random_seed);
0298     }
0299     geo_offset_map.clear();
0300 
0301     for (int layer_i = B0L0_index; layer_i <= B1L1_index; layer_i++)
0302     {
0303         double N_layer_ladder = (layer_i == B0L0_index || layer_i == B0L1_index) ? nLadder_inner : nLadder_outer;
0304 
0305         for (int phi_i = 0; phi_i < N_layer_ladder; phi_i++)
0306         {
0307             geo_offset_map[Form("%i_%i", layer_i, phi_i)] = {
0308                 geo_rand -> Uniform(-random_range_XYZ, random_range_XYZ), 
0309                 geo_rand -> Uniform(-random_range_XYZ, random_range_XYZ),
0310                 geo_rand -> Uniform(-random_range_XYZ, random_range_XYZ)
0311             };
0312         }
0313     }
0314 
0315     geo_offset_correction_fill();
0316 }
0317 
0318 std::vector<double> AvgVtxXY::offset_correction(std::map<std::string,std::vector<double>> input_map)
0319 {
0320     double N_pair = 0;
0321     double sum_x = 0;
0322     double sum_y = 0;
0323     double sum_z = 0;
0324 
0325     for (const auto& pair : input_map)
0326     {
0327         N_pair += 1;
0328         sum_x += pair.second[0];
0329         sum_y += pair.second[1];
0330         sum_z += pair.second[2];
0331     }
0332 
0333     return {sum_x/N_pair, sum_y/N_pair, sum_z/N_pair};
0334 }
0335 
0336 void AvgVtxXY::PreparePairs()
0337 {
0338     cluster_pair_vec.clear();
0339     inner_cluster_vec.clear();
0340     outer_cluster_vec.clear();
0341 
0342     for (int i = 0; i < run_nEvents; i++)
0343     {
0344         tree_in -> GetEntry(i);
0345 
0346         inner_cluster_vec.clear();
0347         outer_cluster_vec.clear();        
0348 
0349         if (MBD_centrality != MBD_centrality) {continue;}
0350         if (is_min_bias != 1) {continue;}
0351         if (MBD_z_vtx != MBD_z_vtx) {continue;}
0352         if (MBD_z_vtx < MBD_vtxZ_cut.first || MBD_z_vtx > MBD_vtxZ_cut.second) {continue;}
0353         if (runnumber != -1 && InttBco_IsToBeRemoved == 1) {continue;}
0354 
0355         if (NClus < INTTNClus_cut.first || NClus > INTTNClus_cut.second) {continue;}
0356 
0357         for (int j = 0; j < NClus; j++)
0358         {
0359             if (ClusAdc->at(j) <= ClusAdc_cut) {continue;}
0360             if (ClusPhiSize->at(j) > ClusPhiSize_cut) {continue;}
0361 
0362             if (geo_offset_map.find(Form("%i_%i", ClusLayer->at(j), ClusLadderPhiId->at(j))) == geo_offset_map.end()) {
0363                 std::cout << "Error : the geo offset map does not have the key : " << Form("%i_%i", ClusLayer->at(j), ClusLadderPhiId->at(j)) << std::endl;
0364                 exit(1);
0365             }
0366 
0367             double post_ClusX = ClusX->at(j) + geo_offset_map[Form("%i_%i", ClusLayer->at(j), ClusLadderPhiId->at(j))][0];
0368             double post_ClusY = ClusY->at(j) + geo_offset_map[Form("%i_%i", ClusLayer->at(j), ClusLadderPhiId->at(j))][1];
0369 
0370             if (ClusLayer->at(j) == B0L0_index || ClusLayer->at(j) == B0L1_index)
0371             {
0372                 inner_cluster_vec.push_back({post_ClusX, post_ClusY});
0373             }
0374             else if (ClusLayer->at(j) == B1L0_index || ClusLayer->at(j) == B1L1_index)
0375             {
0376                 outer_cluster_vec.push_back({post_ClusX, post_ClusY});
0377             }
0378         }
0379 
0380         for (const auto& inner_cluster : inner_cluster_vec)
0381         {
0382             for (const auto& outer_cluster : outer_cluster_vec)
0383             {
0384 
0385                 double Clus_InnerPhiDegree = (inner_cluster.y < 0) ? atan2(inner_cluster.y, inner_cluster.x) * (180./TMath::Pi()) + 360 : atan2(inner_cluster.y, inner_cluster.x) * (180./TMath::Pi()); 
0386                 double Clus_OuterPhiDegree = (outer_cluster.y < 0) ? atan2(outer_cluster.y, outer_cluster.x) * (180./TMath::Pi()) + 360 : atan2(outer_cluster.y, outer_cluster.x) * (180./TMath::Pi()); 
0387 
0388                 // note : to reject the hopeless pairs
0389                 if (fabs(Clus_InnerPhiDegree - Clus_OuterPhiDegree) > hopeless_pair_angle_cut) {continue;}
0390 
0391                 cluster_pair_vec.push_back({inner_cluster, outer_cluster});
0392             }
0393         }
0394 
0395     }
0396 }
0397 
0398 void AvgVtxXY::SetUpHistograms(int iteration, int quadrant)
0399 {
0400     bool insert_check = true;
0401 
0402     // note : inner
0403     insert_check = h2D_map.insert( std::make_pair(
0404         Form("h2_%i_%i_InnerPhi_DeltaPhi", iteration, quadrant), 
0405         new TH2D(
0406             Form("h2_%i_%i_InnerPhi_DeltaPhi", iteration, quadrant),
0407             Form("h2_%i_%i_InnerPhi_DeltaPhi;Inner cluster #phi [radian];#Delta#phi [radian]", iteration, quadrant), 
0408             h2_nbins, Hist2D_angle_xmin, Hist2D_angle_xmax, h2_nbins, Hist2D_DeltaPhi_min, Hist2D_DeltaPhi_max
0409         )
0410     ) ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0411     insert_check = h2D_map.insert( std::make_pair(
0412         Form("h2_%i_%i_InnerPhi_DCA", iteration, quadrant), 
0413         new TH2D(
0414             Form("h2_%i_%i_InnerPhi_DCA", iteration, quadrant),
0415             Form("h2_%i_%i_InnerPhi_DCA;Inner cluster #phi [radian];DCA [cm]", iteration, quadrant), 
0416             h2_nbins, Hist2D_angle_xmin, Hist2D_angle_xmax, h2_nbins, Hist2D_DCA_min, Hist2D_DCA_max
0417         )
0418     ) ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0419 
0420     
0421     // note : outer
0422     insert_check = h2D_map.insert( std::make_pair(
0423         Form("h2_%i_%i_OuterPhi_DeltaPhi", iteration, quadrant), 
0424         new TH2D(
0425             Form("h2_%i_%i_OuterPhi_DeltaPhi", iteration, quadrant),
0426             Form("h2_%i_%i_OuterPhi_DeltaPhi;Outer cluster #phi [radian];#Delta#phi [radian]", iteration, quadrant), 
0427             h2_nbins, Hist2D_angle_xmin, Hist2D_angle_xmax, h2_nbins, Hist2D_DeltaPhi_min, Hist2D_DeltaPhi_max
0428         )
0429     ) ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0430     insert_check = h2D_map.insert( std::make_pair(
0431         Form("h2_%i_%i_OuterPhi_DCA", iteration, quadrant), 
0432         new TH2D(
0433             Form("h2_%i_%i_OuterPhi_DCA", iteration, quadrant),
0434             Form("h2_%i_%i_OuterPhi_DCA;Outer cluster #phi [radian];DCA [cm]", iteration, quadrant), 
0435             h2_nbins, Hist2D_angle_xmin, Hist2D_angle_xmax, h2_nbins, Hist2D_DCA_min, Hist2D_DCA_max
0436         )
0437     ) ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0438     
0439 
0440 
0441     insert_check = h2D_map.insert( std::make_pair(
0442         Form("h2_%i_%i_DeltaPhi_DCA", iteration, quadrant), 
0443         new TH2D(
0444             Form("h2_%i_%i_DeltaPhi_DCA", iteration, quadrant),
0445             Form("h2_%i_%i_DeltaPhi_DCA;#Delta#phi [radian]; DCA [cm]", iteration, quadrant), 
0446             h2_nbins, Hist2D_DeltaPhi_min, Hist2D_DeltaPhi_max, h2_nbins, Hist2D_DCA_min, Hist2D_DCA_max
0447         )
0448     ) ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0449     
0450 
0451     
0452     // note: h1
0453     insert_check = h1D_map.insert( std::make_pair(
0454         Form("h1_%i_%i_DeltaPhi", iteration, quadrant), 
0455         new TH1D(
0456             Form("h1_%i_%i_DeltaPhi", iteration, quadrant),
0457             Form("h1_%i_%i_DeltaPhi;#Delta#phi [radian]; Entries", iteration, quadrant), 
0458             h2_nbins, Hist2D_DeltaPhi_min, Hist2D_DeltaPhi_max
0459         )
0460     ) ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0461 
0462     insert_check = h1D_map.insert( std::make_pair(
0463         Form("h1_%i_%i_DCA", iteration, quadrant), 
0464         new TH1D(
0465             Form("h1_%i_%i_DCA", iteration, quadrant),
0466             Form("h1_%i_%i_DCA;DCA [cm]; Entries", iteration, quadrant), 
0467             h2_nbins, Hist2D_DCA_min, Hist2D_DCA_max
0468         )
0469     ) ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0470 }
0471 
0472 // note : innerphi_delta_phi_fiterror, innerphi_dca_fiterror
0473 std::pair<double,double> AvgVtxXY::RunSingleCandidate(int iteration, int quadrant, std::pair<double,double> assume_center)
0474 {
0475     std::cout<<"iteration : "<<iteration<<" quadrant : "<<quadrant<<", set center: ("<<assume_center.first<<", "<<assume_center.second<<")"<<std::endl;
0476 
0477     SetUpHistograms(iteration, quadrant);
0478 
0479     for (auto &pair : cluster_pair_vec){
0480         
0481         double DCA_sign = calculateAngleBetweenVectors(
0482             pair.second.x, pair.second.y,
0483             pair.first.x, pair.first.y,
0484             assume_center.first, assume_center.second
0485         );
0486 
0487         double Clus_InnerPhi_Offset_radian = atan2(pair.first.y - assume_center.second,  pair.first.x - assume_center.first);
0488         double Clus_OuterPhi_Offset_radian = atan2(pair.second.y - assume_center.second, pair.second.x - assume_center.first);
0489 
0490         double DeltaPhi = Clus_OuterPhi_Offset_radian - Clus_InnerPhi_Offset_radian;
0491 
0492         h2D_map[Form("h2_%i_%i_InnerPhi_DeltaPhi", iteration, quadrant)] -> Fill(Clus_InnerPhi_Offset_radian, DeltaPhi);
0493         h2D_map[Form("h2_%i_%i_InnerPhi_DCA", iteration, quadrant)] -> Fill(Clus_InnerPhi_Offset_radian, DCA_sign);
0494         h2D_map[Form("h2_%i_%i_OuterPhi_DeltaPhi", iteration, quadrant)] -> Fill(Clus_OuterPhi_Offset_radian, DeltaPhi);
0495         h2D_map[Form("h2_%i_%i_OuterPhi_DCA", iteration, quadrant)] -> Fill(Clus_OuterPhi_Offset_radian, DCA_sign);
0496         h2D_map[Form("h2_%i_%i_DeltaPhi_DCA", iteration, quadrant)] -> Fill(DeltaPhi, DCA_sign);
0497 
0498         h1D_map[Form("h1_%i_%i_DeltaPhi", iteration, quadrant)] -> Fill(DeltaPhi);
0499         h1D_map[Form("h1_%i_%i_DCA", iteration, quadrant)] -> Fill(DCA_sign);
0500     }
0501 
0502     bool insert_check = true;
0503 
0504     insert_check = h2D_map.insert( std::make_pair(
0505         Form("h2_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant), 
0506         (TH2D*)h2D_map[Form("h2_%i_%i_InnerPhi_DeltaPhi", iteration, quadrant)]->Clone(Form("h2_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant))
0507     ) ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0508     insert_check = h2D_map.insert( std::make_pair(
0509         Form("h2_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant), 
0510         (TH2D*)h2D_map[Form("h2_%i_%i_InnerPhi_DCA", iteration, quadrant)]->Clone(Form("h2_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant))
0511     ) ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0512 
0513      insert_check = h2D_map.insert( std::make_pair(
0514         Form("h2_%i_%i_OuterPhi_DeltaPhi_bkgrm", iteration, quadrant), 
0515         (TH2D*)h2D_map[Form("h2_%i_%i_OuterPhi_DeltaPhi", iteration, quadrant)]->Clone(Form("h2_%i_%i_OuterPhi_DeltaPhi_bkgrm", iteration, quadrant))
0516     ) ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0517     insert_check = h2D_map.insert( std::make_pair(
0518         Form("h2_%i_%i_OuterPhi_DCA_bkgrm", iteration, quadrant), 
0519         (TH2D*)h2D_map[Form("h2_%i_%i_OuterPhi_DCA", iteration, quadrant)]->Clone(Form("h2_%i_%i_OuterPhi_DCA_bkgrm", iteration, quadrant))
0520     ) ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0521 
0522     TH2D_threshold_advanced_2(h2D_map[Form("h2_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant)], InnerOuterPhi_DeltaPhi_DCA_Threshold);
0523     TH2D_threshold_advanced_2(h2D_map[Form("h2_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant)], InnerOuterPhi_DeltaPhi_DCA_Threshold);
0524     TH2D_threshold_advanced_2(h2D_map[Form("h2_%i_%i_OuterPhi_DeltaPhi_bkgrm", iteration, quadrant)], InnerOuterPhi_DeltaPhi_DCA_Threshold);
0525     TH2D_threshold_advanced_2(h2D_map[Form("h2_%i_%i_OuterPhi_DCA_bkgrm", iteration, quadrant)], InnerOuterPhi_DeltaPhi_DCA_Threshold);
0526     
0527     std::vector<double> h2InnerPhiDeltaPhi_column_content = SumTH2FColumnContent(h2D_map[Form("h2_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant)]);
0528     std::vector<double> h2InnerPhiDCA_column_content = SumTH2FColumnContent(h2D_map[Form("h2_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant)]);
0529 
0530     insert_check = profile_map.insert( std::make_pair(
0531         Form("profile_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant), 
0532         h2D_map[Form("h2_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant)]->ProfileX(Form("profile_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant))
0533     ) ).second; if (insert_check == false) {std::cout<<"Error : profile_map insert failed"<<std::endl; exit(1);}
0534     insert_check = profile_map.insert( std::make_pair(
0535         Form("profile_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant), 
0536         h2D_map[Form("h2_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant)]->ProfileX(Form("profile_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant))
0537     ) ).second; if (insert_check == false) {std::cout<<"Error : profile_map insert failed"<<std::endl; exit(1);}
0538 
0539     insert_check = gr_map.insert( std::make_pair(
0540         Form("gr_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant), 
0541         new TGraph()
0542     ) ).second; if (insert_check == false) {std::cout<<"Error : gr_map insert failed"<<std::endl; exit(1);}
0543     insert_check = gr_map.insert( std::make_pair(
0544         Form("gr_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant), 
0545         new TGraph()
0546     ) ).second; if (insert_check == false) {std::cout<<"Error : gr_map insert failed"<<std::endl; exit(1);}
0547     gr_map[Form("gr_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant)] -> SetName(Form("gr_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant));
0548     gr_map[Form("gr_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant)] -> SetName(Form("gr_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant));
0549 
0550 
0551     for (int i = 0; i < profile_map[Form("profile_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant)]->GetNbinsX(); i++){
0552         if (h2InnerPhiDeltaPhi_column_content[i] < InnerPhi_DeltaPhi_DCA_Fit_Tgr_HitContent_Threshold){continue;} // note : in order to remove some remaining background 
0553 
0554         gr_map[Form("gr_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant)] -> SetPoint(
0555             gr_map[Form("gr_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant)]->GetN(), 
0556             profile_map[Form("profile_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant)]->GetBinCenter(i+1), 
0557             profile_map[Form("profile_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant)]->GetBinContent(i+1)
0558         );
0559     }
0560     for (int i = 0; i < profile_map[Form("profile_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant)]->GetNbinsX(); i++){
0561         if (h2InnerPhiDCA_column_content[i] < InnerPhi_DeltaPhi_DCA_Fit_Tgr_HitContent_Threshold){continue;} // note : in order to remove some remaining background 
0562 
0563         gr_map[Form("gr_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant)] -> SetPoint(
0564             gr_map[Form("gr_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant)]->GetN(), 
0565             profile_map[Form("profile_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant)]->GetBinCenter(i+1), 
0566             profile_map[Form("profile_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant)]->GetBinContent(i+1)
0567         );
0568     }
0569 
0570     // note : do the fit
0571     fit_innerphi_DeltaPhi -> SetParameter(0, 0); 
0572     fit_innerphi_DCA -> SetParameter(0, 0);
0573 
0574     gr_map[Form("gr_%i_%i_InnerPhi_DeltaPhi_bkgrm", iteration, quadrant)] -> Fit(fit_innerphi_DeltaPhi, "NQ", "", -M_PI, M_PI);
0575     gr_map[Form("gr_%i_%i_InnerPhi_DCA_bkgrm", iteration, quadrant)] -> Fit(fit_innerphi_DCA, "NQ", "", -M_PI, M_PI);
0576 
0577     out_iteration = iteration;
0578     out_quadrant = quadrant;
0579     out_assume_center_x = assume_center.first;
0580     out_assume_center_y = assume_center.second;
0581     out_Fit_InnerPhi_DeltaPhi_RChi2 = fit_innerphi_DeltaPhi -> GetChisquare() / double(fit_innerphi_DeltaPhi -> GetNDF()); 
0582     out_Fit_InnerPhi_DeltaPhi_Err0 = fit_innerphi_DeltaPhi -> GetParError(0);
0583     out_Fit_InnerPhi_DeltaPhi_LineY = fit_innerphi_DeltaPhi -> GetParameter(0);
0584     out_Fit_InnerPhi_DCA_RChi2 = fit_innerphi_DCA -> GetChisquare() / double(fit_innerphi_DCA -> GetNDF());
0585     out_Fit_InnerPhi_DCA_Err0 = fit_innerphi_DCA -> GetParError(0);
0586     out_Fit_InnerPhi_DCA_LineY = fit_innerphi_DCA -> GetParameter(0);
0587 
0588     tree_quadrant_detail -> Fill();
0589 
0590     std::cout<<"iteration : "<<iteration<<" quadrant : "<<quadrant<<", FitErrors: ("<<fit_innerphi_DeltaPhi -> GetParError(0)<<", "<<fit_innerphi_DCA -> GetParError(0)<<")"<<std::endl;
0591 
0592     return{
0593         fit_innerphi_DeltaPhi -> GetParError(0),
0594         fit_innerphi_DCA -> GetParError(0)
0595     };
0596 
0597 }
0598 
0599 void AvgVtxXY::FindVertexQuadrant(int nIteration, double search_Window_half, std::pair<double,double> set_center)
0600 {
0601     double dynamic_search_window_half = search_Window_half;
0602     std::pair<double, double> dynamic_center = set_center;
0603 
0604     int small_index;
0605     std::pair<double,double> small_pair;
0606 
0607     std::vector<std::pair<double,double>> vtx_vec = Get4vtx(dynamic_center,dynamic_search_window_half);
0608 
0609     bool insert_check = true;
0610 
0611     insert_check = gr_map.insert( std::make_pair(
0612         Form("History_InnerPhi_DeltaPhi_fitErrorAll"), 
0613         new TGraph()
0614     ) ).second; if (insert_check == false) {std::cout<<"Error : gr_map insert failed"<<std::endl; exit(1);}
0615     insert_check = gr_map.insert( std::make_pair(
0616         Form("History_InnerPhi_DeltaPhi_fitErrorWinner"), 
0617         new TGraph()
0618     ) ).second; if (insert_check == false) {std::cout<<"Error : gr_map insert failed"<<std::endl; exit(1);}
0619 
0620     insert_check = gr_map.insert( std::make_pair(
0621         Form("History_InnerPhi_DCA_fitErrorAll"), 
0622         new TGraph()
0623     ) ).second; if (insert_check == false) {std::cout<<"Error : gr_map insert failed"<<std::endl; exit(1);}
0624     insert_check = gr_map.insert( std::make_pair(
0625         Form("History_InnerPhi_DCA_fitErrorWinner"), 
0626         new TGraph()
0627     ) ).second; if (insert_check == false) {std::cout<<"Error : gr_map insert failed"<<std::endl; exit(1);}
0628 
0629     insert_check = gr_map.insert( std::make_pair(
0630         Form("History_Winner_index"), 
0631         new TGraph()
0632     ) ).second; if (insert_check == false) {std::cout<<"Error : gr_map insert failed"<<std::endl; exit(1);}
0633 
0634     gr_map[Form("History_InnerPhi_DeltaPhi_fitErrorAll")]->SetName(Form("History_InnerPhi_DeltaPhi_fitErrorAll"));
0635     gr_map[Form("History_InnerPhi_DeltaPhi_fitErrorWinner")]->SetName(Form("History_InnerPhi_DeltaPhi_fitErrorWinner"));
0636     gr_map[Form("History_InnerPhi_DCA_fitErrorAll")]->SetName(Form("History_InnerPhi_DCA_fitErrorAll"));
0637     gr_map[Form("History_InnerPhi_DCA_fitErrorWinner")]->SetName(Form("History_InnerPhi_DCA_fitErrorWinner"));
0638     gr_map[Form("History_Winner_index")]->SetName(Form("History_Winner_index"));
0639 
0640 
0641     for (int trial_i = 0; trial_i < nIteration; trial_i++)
0642     {   
0643         for (int quadrant_i = 0; quadrant_i < 4; quadrant_i++)
0644         {
0645             // note : innerphi_delta_phi_fiterror, innerphi_dca_fiterror
0646             std::pair<double,double> current_pair = RunSingleCandidate(trial_i, quadrant_i, vtx_vec[quadrant_i]);
0647 
0648             gr_map[Form("History_InnerPhi_DeltaPhi_fitErrorAll")] -> SetPoint(
0649                 gr_map[Form("History_InnerPhi_DeltaPhi_fitErrorAll")]->GetN(), 
0650                 trial_i, 
0651                 current_pair.first
0652             );
0653 
0654             gr_map[Form("History_InnerPhi_DCA_fitErrorAll")] -> SetPoint(
0655                 gr_map[Form("History_InnerPhi_DCA_fitErrorAll")]->GetN(), 
0656                 trial_i, 
0657                 current_pair.second
0658             );
0659 
0660             if (quadrant_i == 0)
0661             {
0662                 small_index = quadrant_i;
0663                 small_pair = current_pair;
0664             }
0665             else 
0666             {
0667                 if (current_pair.first < small_pair.first && current_pair.second < small_pair.second)
0668                 // if (current_pair.second < small_pair.second) // note : only consider the DCA fit error
0669                 {
0670                     small_index = quadrant_i;
0671                     small_pair = current_pair;
0672                 }
0673             }
0674         }
0675 
0676         gr_map[Form("History_Winner_index")] -> SetPoint(
0677             gr_map[Form("History_Winner_index")]->GetN(), 
0678             trial_i, 
0679             small_index
0680         );
0681 
0682         gr_map[Form("History_InnerPhi_DeltaPhi_fitErrorWinner")] -> SetPoint(
0683             gr_map[Form("History_InnerPhi_DeltaPhi_fitErrorWinner")]->GetN(), 
0684             trial_i, 
0685             small_pair.first
0686         );
0687 
0688         gr_map[Form("History_InnerPhi_DCA_fitErrorWinner")] -> SetPoint(
0689             gr_map[Form("History_InnerPhi_DCA_fitErrorWinner")]->GetN(), 
0690             trial_i, 
0691             small_pair.second
0692         );
0693 
0694         std::cout<<"In FindVertexQuadrant, iteration : "<<trial_i<<" small_index : "<<small_index<<" small_pair : ("<<vtx_vec[small_index].first<<", "<<vtx_vec[small_index].second<<")"<<std::endl;
0695 
0696         // note : generating the new 4 vertex for the next comparison
0697         // note : start to shrink the square
0698         if (trial_i != nIteration - 1)
0699         {
0700             dynamic_center = {(vtx_vec[small_index].first + dynamic_center.first)/2., (vtx_vec[small_index].second + dynamic_center.second)/2.};
0701             // cout<<"test : "<<origin.first<<" "<<origin.second<<" length: "<<length<<endl;
0702             dynamic_search_window_half /= 2.;
0703             vtx_vec = Get4vtx(dynamic_center,dynamic_search_window_half); 
0704         }
0705 
0706     }
0707 
0708     final_InnerPhi_DeltaPhi = (TH2D*)h2D_map[Form("h2_%i_%i_InnerPhi_DeltaPhi", nIteration-1, small_index)]->Clone("final_InnerPhi_DeltaPhi");
0709     final_InnerPhi_DCA = (TH2D*)h2D_map[Form("h2_%i_%i_InnerPhi_DCA", nIteration-1, small_index)]->Clone("final_InnerPhi_DCA");
0710     final_InnerPhi_DeltaPhi_bkgrm = (TH2D*)h2D_map[Form("h2_%i_%i_InnerPhi_DeltaPhi_bkgrm", nIteration-1, small_index)]->Clone("final_InnerPhi_DeltaPhi_bkgrm");
0711     final_InnerPhi_DCA_bkgrm = (TH2D*)h2D_map[Form("h2_%i_%i_InnerPhi_DCA_bkgrm", nIteration-1, small_index)]->Clone("final_InnerPhi_DCA_bkgrm");
0712 
0713     final_OuterPhi_DeltaPhi = (TH2D*)h2D_map[Form("h2_%i_%i_OuterPhi_DeltaPhi", nIteration-1, small_index)]->Clone("final_OuterPhi_DeltaPhi");
0714     final_OuterPhi_DCA = (TH2D*)h2D_map[Form("h2_%i_%i_OuterPhi_DCA", nIteration-1, small_index)]->Clone("final_OuterPhi_DCA");
0715     final_OuterPhi_DeltaPhi_bkgrm = (TH2D*)h2D_map[Form("h2_%i_%i_OuterPhi_DeltaPhi_bkgrm", nIteration-1, small_index)]->Clone("final_OuterPhi_DeltaPhi_bkgrm");
0716     final_OuterPhi_DCA_bkgrm = (TH2D*)h2D_map[Form("h2_%i_%i_OuterPhi_DCA_bkgrm", nIteration-1, small_index)]->Clone("final_OuterPhi_DCA_bkgrm");
0717 
0718     final_DeltaPhi_DCA = (TH2D*)h2D_map[Form("h2_%i_%i_DeltaPhi_DCA", nIteration-1, small_index)]->Clone("final_DeltaPhi_DCA");
0719     final_DeltaPhi_1D = (TH1D*)h1D_map[Form("h1_%i_%i_DeltaPhi", nIteration-1, small_index)]->Clone("final_DeltaPhi");
0720     final_DCA_1D = (TH1D*)h1D_map[Form("h1_%i_%i_DCA", nIteration-1, small_index)]->Clone("final_DCA");
0721 
0722 
0723     out_vtxX_quadrant = (vtx_vec[small_index].first + dynamic_center.first)/2.;
0724     out_vtxXE_quadrant = fabs(vtx_vec[small_index].first - dynamic_center.first)/2.;
0725     out_vtxY_quadrant = (vtx_vec[small_index].second + dynamic_center.second)/2.;;
0726     out_vtxYE_quadrant = fabs(vtx_vec[small_index].second - dynamic_center.second)/2.;    
0727 
0728 }
0729 
0730 void AvgVtxXY::FindVertexLineFill(std::pair<double,double> set_center, int Nbins, double search_Window_half, double segmentation)
0731 {
0732     bool insert_check = true;
0733 
0734     insert_check = h2D_map.insert(
0735         std::make_pair(
0736             "h2_LineFilled_XY",
0737             new TH2D(
0738                 "h2_LineFilled_XY",
0739                 "h2_LineFilled_XY;X [cm];Y [cm]",
0740                 Nbins, set_center.first - search_Window_half, set_center.first + search_Window_half,
0741                 Nbins, set_center.second - search_Window_half, set_center.second + search_Window_half
0742             )
0743         )
0744     ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0745 
0746     // insert_check = h2D_map.insert(
0747     //     std::make_pair(
0748     //         "h2_LineFilled_XY_bkgrm",
0749     //         new TH2D(
0750     //             "h2_LineFilled_XY_bkgrm",
0751     //             "h2_LineFilled_XY_bkgrm;X [cm];Y [cm]",
0752     //             Nbins, set_center.first - search_Window_half, set_center.first + search_Window_half,
0753     //             Nbins, set_center.second - search_Window_half, set_center.second + search_Window_half
0754     //         )
0755     //     )
0756     // ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0757 
0758     for (int i = 0; i < cluster_pair_vec.size(); i++){
0759         double Clus_InnerPhi_Offset = (cluster_pair_vec[i].first.y - set_center.second < 0) ? atan2(cluster_pair_vec[i].first.y - set_center.second, cluster_pair_vec[i].first.x - set_center.first) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].first.y - set_center.second, cluster_pair_vec[i].first.x - set_center.first) * (180./TMath::Pi());
0760         double Clus_OuterPhi_Offset = (cluster_pair_vec[i].second.y - set_center.second < 0) ? atan2(cluster_pair_vec[i].second.y - set_center.second, cluster_pair_vec[i].second.x - set_center.first) * (180./TMath::Pi()) + 360 : atan2(cluster_pair_vec[i].second.y - set_center.second, cluster_pair_vec[i].second.x - set_center.first) * (180./TMath::Pi());
0761 
0762         std::vector<double> DCA_info_vec = calculateDistanceAndClosestPoint(
0763             cluster_pair_vec[i].first.x, cluster_pair_vec[i].first.y,
0764             cluster_pair_vec[i].second.x, cluster_pair_vec[i].second.y,
0765             set_center.first, set_center.second
0766         );
0767 
0768         if (fabs(Clus_InnerPhi_Offset - Clus_OuterPhi_Offset) < rough_DeltaPhi_cut_degree)
0769         {
0770             TH2DSampleLineFill(h2D_map["h2_LineFilled_XY"], segmentation, {cluster_pair_vec[i].first.x, cluster_pair_vec[i].first.y}, {DCA_info_vec[1], DCA_info_vec[2]});
0771             // std::cout<<"filled checked, cluster pos : "<<cluster_pair_vec[i].first.x<<" "<<cluster_pair_vec[i].first.y<<" "<<cluster_pair_vec[i].second.x<<" "<<cluster_pair_vec[i].second.y<<std::endl;
0772         }
0773     }
0774     
0775     insert_check = h2D_map.insert(
0776         std::make_pair(
0777             "h2_LineFilled_XY_bkgrm",
0778             (TH2D*) h2D_map["h2_LineFilled_XY"]->Clone("h2_LineFilled_XY_bkgrm")
0779         )
0780     ).second; if (insert_check == false) {std::cout<<"Error : h2D_map insert failed"<<std::endl; exit(1);}
0781 
0782     TH2D_threshold_advanced_2(h2D_map["h2_LineFilled_XY_bkgrm"], LineFill_Hist_Threshold);
0783 
0784     out_vtxX_LineFill       = h2D_map["h2_LineFilled_XY_bkgrm"] -> GetMean(1);     // + xy_hist_bkgrm -> GetXaxis() -> GetBinWidth(1) / 2.; // note : the TH2F calculate the GetMean based on the bin center, no need to apply additional offset
0785     out_vtxXE_LineFill      = h2D_map["h2_LineFilled_XY_bkgrm"] -> GetMeanError(1);
0786     out_vtxXStdDev_LineFill = h2D_map["h2_LineFilled_XY_bkgrm"] -> GetStdDev(1);
0787     out_vtxY_LineFill       = h2D_map["h2_LineFilled_XY_bkgrm"] -> GetMean(2);     // + xy_hist_bkgrm -> GetYaxis() -> GetBinWidth(1) / 2.; // note : the TH2F calculate the GetMean based on the bin center, no need to apply additional offset
0788     out_vtxYE_LineFill      = h2D_map["h2_LineFilled_XY_bkgrm"] -> GetMeanError(2);
0789     out_vtxYStdDev_LineFill = h2D_map["h2_LineFilled_XY_bkgrm"] -> GetStdDev(2);
0790 
0791     return;
0792 }
0793 
0794 void AvgVtxXY::TH2DSampleLineFill(TH2D * hist_in, double segmentation, std::pair<double,double> inner_clu, std::pair<double,double> outer_clu)
0795 {
0796     double x_min = hist_in -> GetXaxis() -> GetXmin();
0797     double x_max = hist_in -> GetXaxis() -> GetXmax();
0798     double y_min = hist_in -> GetYaxis() -> GetXmin();
0799     double y_max = hist_in -> GetYaxis() -> GetXmax();
0800 
0801     double seg_x, seg_y;
0802     double angle;
0803     int n_seg = 0;
0804 
0805     while (true)
0806     {
0807         angle = atan2(inner_clu.second-outer_clu.second, inner_clu.first-outer_clu.first);
0808         seg_x = (n_seg * segmentation) * cos(angle) + outer_clu.first; // note : atan2(y,x), point.first is the radius
0809         seg_y = (n_seg * segmentation) * sin(angle) + outer_clu.second;
0810         
0811         if ( (seg_x > x_min && seg_x < x_max && seg_y > y_min && seg_y < y_max) != true ) {break;}
0812         hist_in -> Fill(seg_x, seg_y);
0813         n_seg += 1;
0814     }
0815 
0816     n_seg = 1;
0817     while (true)
0818     {
0819         angle = atan2(inner_clu.second-outer_clu.second, inner_clu.first-outer_clu.first);
0820         seg_x = (-1 * n_seg * segmentation) * cos(angle) + outer_clu.first; // note : atan2(y,x), point.first is the radius
0821         seg_y = (-1 * n_seg * segmentation) * sin(angle) + outer_clu.second;
0822         
0823         if ( (seg_x > x_min && seg_x < x_max && seg_y > y_min && seg_y < y_max) != true ) {break;}
0824         hist_in -> Fill(seg_x, seg_y);
0825         n_seg += 1;
0826     }
0827 }
0828 
0829 std::vector<double> AvgVtxXY::calculateDistanceAndClosestPoint(double x1, double y1, double x2, double y2, double target_x, double target_y) 
0830 {
0831     
0832     if (x1 != x2)
0833     {
0834         // Calculate the slope and intercept of the line passing through (x1, y1) and (x2, y2)
0835         double a = (y2 - y1) / (x2 - x1);
0836         double b = y1 - a * x1;
0837 
0838         // cout<<"slope : y="<<a<<"x+"<<b<<endl;
0839         
0840         // Calculate the closest distance from (target_x, target_y) to the line y = ax + b
0841         double closest_distance = std::abs(a * target_x - target_y + b) / std::sqrt(a * a + 1);
0842 
0843         // Calculate the coordinates of the closest point (Xc, Yc) on the line y = ax + b
0844         double Xc = (target_x + a * target_y - a * b) / (a * a + 1);
0845         double Yc = a * Xc + b;
0846 
0847         return { closest_distance, Xc, Yc };
0848     }
0849     else 
0850     {
0851         double closest_distance = std::abs(x1 - target_x);
0852         double Xc = x1;
0853         double Yc = target_y;
0854 
0855         return { closest_distance, Xc, Yc };
0856     }
0857     
0858     
0859 }
0860 
0861 void AvgVtxXY::EndRun()
0862 {
0863     std::cout<<std::endl;
0864     std::cout<<"In runnumber : "<<runnumber<<std::endl;
0865     std::cout<<"In job : "<<process_id<<std::endl;
0866     std::cout<<"w/ Nevent: "<<run_nEvents<<std::endl;
0867     std::cout<<"w/ NCluster pair: "<<cluster_pair_vec.size()<<std::endl;
0868     std::cout<<"Quadrant method: ("<<Form("%.5f +- %.5f [cm], %.5f +- %.5f [cm]", out_vtxX_quadrant, out_vtxXE_quadrant, out_vtxY_quadrant, out_vtxYE_quadrant)<<")"<<std::endl;
0869     std::cout<<"LineFill method: ("<<Form("%.5f +- %.5f [cm], %.5f +- %.5f [cm]", out_vtxX_LineFill, out_vtxXE_LineFill, out_vtxY_LineFill, out_vtxYE_LineFill)<<")"<<std::endl;
0870     std::cout<<"Final : "<<Form("(%.5f [cm], %.5f [cm])", (out_vtxX_quadrant + out_vtxX_LineFill)/2., (out_vtxY_quadrant + out_vtxY_LineFill)/2.)<<std::endl;
0871 
0872     file_out -> cd();
0873     tree_vtxXY -> Fill();
0874     tree_vtxXY -> Write();
0875     if (HaveGeoOffsetTag) {tree_geooffset -> Write();}
0876     tree_quadrant_detail -> Write();
0877     tree_parameters -> Write();
0878 
0879     for (const auto& pair : h2D_map)
0880     {
0881         pair.second -> Write();
0882     }
0883 
0884     for (const auto& pair : h1D_map)
0885     {
0886         pair.second -> Write();
0887     }
0888 
0889     for (const auto& pair : gr_map)
0890     {
0891         pair.second -> Write();
0892     }
0893 
0894     for (const auto& pair : profile_map)
0895     {
0896         pair.second -> Write();
0897     }
0898 
0899     final_InnerPhi_DeltaPhi -> Write();
0900     final_InnerPhi_DCA -> Write();
0901     final_InnerPhi_DeltaPhi_bkgrm -> Write();
0902     final_InnerPhi_DCA_bkgrm -> Write();
0903     final_OuterPhi_DeltaPhi -> Write();
0904     final_OuterPhi_DCA -> Write();
0905     final_OuterPhi_DeltaPhi_bkgrm -> Write();
0906     final_OuterPhi_DCA_bkgrm -> Write();
0907     final_DeltaPhi_DCA -> Write();
0908     final_DeltaPhi_1D -> Write();
0909     final_DCA_1D -> Write();
0910 
0911     file_out -> Close();
0912 
0913 }
0914 
0915 double AvgVtxXY::calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY) {
0916     // Calculate the vectors vector_1 (point_1 to point_2) and vector_2 (point_1 to target)
0917     double vector1X = x2 - x1;
0918     double vector1Y = y2 - y1;
0919 
0920     double vector2X = targetX - x1;
0921     double vector2Y = targetY - y1;
0922 
0923     // Calculate the cross product of vector_1 and vector_2 (z-component)
0924     double crossProduct = vector1X * vector2Y - vector1Y * vector2X;
0925     
0926     // cout<<" crossProduct : "<<crossProduct<<endl;
0927 
0928     // Calculate the magnitudes of vector_1 and vector_2
0929     double magnitude1 = std::sqrt(vector1X * vector1X + vector1Y * vector1Y);
0930     double magnitude2 = std::sqrt(vector2X * vector2X + vector2Y * vector2Y);
0931 
0932     // Calculate the angle in radians using the inverse tangent of the cross product and dot product
0933     double dotProduct = vector1X * vector2X + vector1Y * vector2Y;
0934 
0935     double angleInRadians = std::atan2(std::abs(crossProduct), dotProduct);
0936     // Convert the angle from radians to degrees and return it
0937     double angleInDegrees = angleInRadians * 180.0 / M_PI;
0938     
0939     double angleInRadians_new = std::asin( crossProduct/(magnitude1*magnitude2) );
0940     double angleInDegrees_new = angleInRadians_new * 180.0 / M_PI;
0941     
0942     // cout<<"angle : "<<angleInDegrees_new<<endl;
0943 
0944     double DCA_distance = sin(angleInRadians_new) * magnitude2;
0945 
0946     return DCA_distance;
0947 }
0948 
0949 void AvgVtxXY::TH2D_threshold_advanced_2(TH2D * hist, double threshold)
0950 {
0951     // note : this function is to remove the background of the 2D histogram
0952     // note : but the threshold is given by average of the contents of the top "chosen_bin" bins and timing the threshold
0953     // todo : here
0954     double max_cut = 0;
0955     int chosen_bin = TH2D_threshold_advanced_2_CheckNBins;
0956 
0957     std::vector<float> all_bin_content_vec; all_bin_content_vec.clear();
0958     for (int xi = 0; xi < hist -> GetNbinsX(); xi++){
0959         for(int yi = 0; yi < hist -> GetNbinsY(); yi++){
0960             all_bin_content_vec.push_back(hist -> GetBinContent(xi + 1, yi +1));
0961         }
0962     }
0963     std::vector<unsigned long> ind(all_bin_content_vec.size(),0);
0964     TMath::Sort(all_bin_content_vec.size(), &all_bin_content_vec[0], &ind[0]);
0965     for (int i = 0; i < chosen_bin; i++) {max_cut += all_bin_content_vec[ind[i]]; /*cout<<"test : "<<all_bin_content_vec[ind[i]]<<endl;*/}
0966 
0967     max_cut = (max_cut / double(chosen_bin)) * threshold;
0968 
0969     for (int xi = 0; xi < hist -> GetNbinsX(); xi++){
0970         for(int yi = 0; yi < hist -> GetNbinsY(); yi++){
0971             if (hist -> GetBinContent(xi + 1, yi +1) < max_cut){ hist -> SetBinContent(xi + 1, yi +1, 0); }
0972         }
0973     }
0974 }
0975 
0976 std::vector<double> AvgVtxXY::SumTH2FColumnContent(TH2D * hist_in)
0977 {
0978     std::vector<double> sum_vec; sum_vec.clear();
0979     for (int i = 0; i < hist_in -> GetNbinsX(); i++){
0980         double sum = 0;
0981         for (int j = 0; j < hist_in -> GetNbinsY(); j++){
0982             sum += hist_in -> GetBinContent(i+1, j+1);
0983         }
0984         sum_vec.push_back(sum);
0985     }
0986     return sum_vec;
0987 }
0988 
0989 std::vector<std::pair<double,double>> AvgVtxXY::Get4vtx(std::pair<double,double> origin, double length)
0990 {
0991     std::vector<std::pair<double,double>> unit_vtx = {{1,1},{-1,1},{-1,-1},{1,-1}};
0992     std::vector<std::pair<double,double>> vec_out; vec_out.clear();
0993 
0994     for (std::pair i1 : unit_vtx)
0995     {
0996         vec_out.push_back({i1.first * length + origin.first, i1.second * length + origin.second});
0997     }
0998 
0999     return vec_out;
1000 }
[ 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!