Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-09 08:12:15 

0001 #include "RestDist.h"
0002 
0003 RestDist::RestDist(
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 
0011     std::string output_file_name_suffix_in,
0012     std::pair<double, double> vertexXYIncm_in,
0013 
0014     bool Apply_cut_in,
0015     bool ApplyVtxZReWeighting_in,
0016     std::pair<bool, int> ApplyEvtBcoFullDiffCut_in,
0017     std::pair<bool, std::pair<double,double>> RequireVtxZRange_in,
0018     std::pair<bool, std::pair<double,double>> isClusQA_in, // note : true/false, {ADC, phi_size}
0019     bool isRotated_in
0020 ):
0021     process_id(process_id_in),
0022     runnumber(runnumber_in),
0023     run_nEvents(run_nEvents_in),
0024     input_directory(input_directory_in),
0025     input_file_name(input_file_name_in),
0026     output_directory(output_directory_in),
0027     output_file_name_suffix(output_file_name_suffix_in),
0028     vertexXYIncm(vertexXYIncm_in),
0029     Apply_cut(Apply_cut_in),
0030     ApplyVtxZReWeighting(ApplyVtxZReWeighting_in),
0031     ApplyEvtBcoFullDiffCut(ApplyEvtBcoFullDiffCut_in),
0032     RequireVtxZRange(RequireVtxZRange_in),
0033     isClusQA(isClusQA_in),
0034     isRotated(isRotated_in)
0035 {
0036     PrepareOutPutFileName();
0037     PrepareInputFile();
0038 
0039     run_nEvents = (run_nEvents == -1) ? tree_in->GetEntries() : run_nEvents;
0040     run_nEvents = (run_nEvents > tree_in->GetEntries()) ? tree_in->GetEntries() : run_nEvents;
0041     nCentrality_bin = centrality_edges.size() - 1;
0042 
0043     PrepareOutputFile();
0044     PrepareHist();
0045 
0046     temp_all_pair_vec.clear();
0047     temp_all_pair_deltaR_vec.clear();
0048     BestPair_h2D_map.clear();
0049     BestPair_gr_map.clear();
0050 
0051 }
0052 
0053 void RestDist::PrepareOutPutFileName()
0054 {
0055     std::string job_index = std::to_string( process_id );
0056     int job_index_len = 5;
0057     job_index.insert(0, job_index_len - job_index.size(), '0');
0058 
0059     std::string runnumber_str = std::to_string( runnumber );
0060     if (runnumber != -1){
0061         int runnumber_str_len = 8;
0062         runnumber_str.insert(0, runnumber_str_len - runnumber_str.size(), '0');
0063     }
0064 
0065     if (output_file_name_suffix.size() > 0 && output_file_name_suffix[0] != '_') {
0066         output_file_name_suffix = "_" + output_file_name_suffix;
0067     }
0068 
0069     output_filename = "RestDist";
0070     output_filename = (runnumber != -1) ? "Data_" + output_filename : "MC_" + output_filename;
0071 
0072     output_filename += (Apply_cut) ? "_vtxZQA" : "_NovtxZQA";
0073     output_filename += (ApplyVtxZReWeighting) ? "_VtxZReWeighting" : "";
0074     output_filename += (runnumber != -1 && ApplyEvtBcoFullDiffCut.first) ? Form("_EvtBcoFullDiffCut%d", ApplyEvtBcoFullDiffCut.second) : "";
0075     output_filename += (isRotated) ? "_Rotated" : "";
0076 
0077     if (RequireVtxZRange.first){
0078         std::string vtxZRange;
0079 
0080         double vtxZ_range_l = RequireVtxZRange.second.first;
0081         double vtxZ_range_r = RequireVtxZRange.second.second;
0082 
0083         vtxZRange += "_vtxZRange";
0084         vtxZRange += (vtxZ_range_l < 0) ? Form("M%.0fp%d",fabs(vtxZ_range_l),int(fabs(vtxZ_range_l)*10)%10) : Form("%.0fp%d",vtxZ_range_l, int(vtxZ_range_l*10)%10);
0085         vtxZRange += "to";
0086         vtxZRange += (vtxZ_range_r < 0) ? Form("M%.0fp%d",fabs(vtxZ_range_r),int(fabs(vtxZ_range_r)*10)%10) : Form("%.0fp%d",vtxZ_range_r, int(vtxZ_range_r*10)%10);
0087 
0088         output_filename += vtxZRange;
0089     }
0090 
0091     output_filename += (isClusQA.first) ? Form("_ClusQAAdc%.0fPhiSize%.0f",isClusQA.second.first,isClusQA.second.second) : "";
0092 
0093     output_filename += output_file_name_suffix;
0094     output_filename += (runnumber != -1) ? Form("_%s_%s.root",runnumber_str.c_str(),job_index.c_str()) : Form("_%s.root",job_index.c_str());
0095 }
0096 
0097 std::map<std::string, int> RestDist::GetInputTreeBranches(TTree * m_tree_in)
0098 {
0099     std::map<std::string, int> branch_map;
0100     TObjArray * branch_list = m_tree_in -> GetListOfBranches();
0101     for (int i = 0; i < branch_list -> GetEntries(); i++)
0102     {
0103         TBranch * branch = dynamic_cast<TBranch*>(branch_list->At(i));
0104         branch_map[branch -> GetName()] = 1;
0105     }
0106     return branch_map;
0107 }
0108 
0109 void RestDist::PrepareInputFile()
0110 {
0111     file_in = TFile::Open(Form("%s/%s", input_directory.c_str(), input_file_name.c_str()));
0112     tree_in = (TTree*)file_in->Get(tree_name.c_str());
0113 
0114     std::map<std::string, int> branch_map = GetInputTreeBranches(tree_in);
0115 
0116     tree_in -> SetBranchStatus("*", 0);
0117 
0118     tree_in -> SetBranchStatus("is_min_bias", 1);
0119     tree_in -> SetBranchStatus("MBD_centrality", 1);
0120     tree_in -> SetBranchStatus("MBD_z_vtx", 1);
0121     tree_in -> SetBranchStatus("MBD_charge_sum", 1);
0122     tree_in -> SetBranchStatus("MBD_charge_asymm", 1);
0123     tree_in -> SetBranchStatus("MBD_south_charge_sum", 1);
0124     tree_in -> SetBranchStatus("MBD_north_charge_sum", 1);
0125 
0126     tree_in -> SetBranchStatus("INTTvtxZ", 1);
0127     tree_in -> SetBranchStatus("INTTvtxZError", 1);
0128     tree_in -> SetBranchStatus("NgroupTrapezoidal", 1);
0129     tree_in -> SetBranchStatus("NgroupCoarse", 1);
0130     tree_in -> SetBranchStatus("TrapezoidalFitWidth", 1);
0131     tree_in -> SetBranchStatus("TrapezoidalFWHM", 1);
0132 
0133     tree_in -> SetBranchStatus("NClus", 1);
0134     tree_in -> SetBranchStatus("NClus_Layer1", 1);
0135     
0136     // note : for data
0137     if (branch_map.find("MBDNSg2") != branch_map.end()) {
0138         tree_in -> SetBranchStatus("MBDNSg2", 1);
0139         m_withTrig = true;
0140     }
0141     if (branch_map.find("MBDNSg2_vtxZ10cm") != branch_map.end()) {tree_in -> SetBranchStatus("MBDNSg2_vtxZ10cm", 1);}
0142     if (branch_map.find("MBDNSg2_vtxZ30cm") != branch_map.end()) {tree_in -> SetBranchStatus("MBDNSg2_vtxZ30cm", 1);}
0143     
0144     if (branch_map.find("InttBcoFullDiff_next") != branch_map.end()) {tree_in -> SetBranchStatus("InttBcoFullDiff_next", 1); }
0145     
0146     // note :for MC
0147     if (branch_map.find("NTruthVtx") != branch_map.end()) {tree_in -> SetBranchStatus("NTruthVtx", 1);}
0148     if (branch_map.find("TruthPV_trig_x") != branch_map.end()) {tree_in -> SetBranchStatus("TruthPV_trig_x", 1);}
0149     if (branch_map.find("TruthPV_trig_y") != branch_map.end()) {tree_in -> SetBranchStatus("TruthPV_trig_y", 1);}
0150     if (branch_map.find("TruthPV_trig_z") != branch_map.end()) {tree_in -> SetBranchStatus("TruthPV_trig_z", 1);}
0151 
0152 
0153     tree_in -> SetBranchStatus("ClusX", 1);
0154     tree_in -> SetBranchStatus("ClusY", 1);
0155     tree_in -> SetBranchStatus("ClusZ", 1);
0156 
0157     tree_in -> SetBranchStatus("ClusAdc", 1);
0158     tree_in -> SetBranchStatus("ClusPhiSize", 1);
0159     tree_in -> SetBranchStatus("ClusZSize", 1);
0160 
0161     tree_in -> SetBranchStatus("ClusLayer", 1);
0162     tree_in -> SetBranchStatus("ClusLadderZId", 1);
0163     tree_in -> SetBranchStatus("ClusLadderPhiId", 1);
0164 
0165     tree_in -> SetBranchStatus("ClusEta_INTTz", 1);
0166     tree_in -> SetBranchStatus("ClusEta_MBDz", 1);
0167     tree_in -> SetBranchStatus("ClusPhi_AvgPV", 1);
0168 
0169     // note : for MC
0170     if (branch_map.find("ClusEta_TrueXYZ") != branch_map.end()) {tree_in -> SetBranchStatus("ClusEta_TrueXYZ", 1);}
0171     if (branch_map.find("ClusPhi_TrueXY") != branch_map.end()) {tree_in -> SetBranchStatus("ClusPhi_TrueXY", 1);}
0172 
0173     // Division : ---SetBranchAddress-----------------------------------------------------------------------------------------------
0174 
0175 
0176     ClusX = 0;
0177     ClusY = 0;
0178     ClusZ = 0;
0179     ClusAdc = 0;
0180     ClusPhiSize = 0;
0181     ClusZSize = 0;
0182     ClusLayer = 0;
0183     ClusLadderZId = 0;
0184     ClusLadderPhiId = 0;
0185     ClusEta_INTTz = 0;
0186     ClusEta_MBDz = 0;
0187     ClusPhi_AvgPV = 0;
0188     ClusEta_TrueXYZ = 0;
0189     ClusPhi_TrueXY = 0;
0190 
0191     tree_in -> SetBranchAddress("is_min_bias", &is_min_bias);
0192     tree_in -> SetBranchAddress("MBD_centrality", &MBD_centrality);
0193     tree_in -> SetBranchAddress("MBD_z_vtx", &MBD_z_vtx);
0194     tree_in -> SetBranchAddress("MBD_charge_sum", &MBD_charge_sum);
0195     tree_in -> SetBranchAddress("MBD_charge_asymm", &MBD_charge_asymm);
0196     tree_in -> SetBranchAddress("MBD_south_charge_sum", &MBD_south_charge_sum);
0197     tree_in -> SetBranchAddress("MBD_north_charge_sum", &MBD_north_charge_sum);
0198 
0199     tree_in -> SetBranchAddress("INTTvtxZ", &INTTvtxZ);
0200     tree_in -> SetBranchAddress("INTTvtxZError", &INTTvtxZError);
0201     tree_in -> SetBranchAddress("NgroupTrapezoidal", &NgroupTrapezoidal);
0202     tree_in -> SetBranchAddress("NgroupCoarse", &NgroupCoarse);
0203     tree_in -> SetBranchAddress("TrapezoidalFitWidth", &TrapezoidalFitWidth);
0204     tree_in -> SetBranchAddress("TrapezoidalFWHM", &TrapezoidalFWHM);
0205 
0206     tree_in -> SetBranchAddress("NClus", &NClus);
0207     tree_in -> SetBranchAddress("NClus_Layer1", &NClus_Layer1);
0208 
0209     tree_in -> SetBranchAddress("ClusX", &ClusX);
0210     tree_in -> SetBranchAddress("ClusY", &ClusY);
0211     tree_in -> SetBranchAddress("ClusZ", &ClusZ);
0212 
0213     tree_in -> SetBranchAddress("ClusAdc", &ClusAdc);
0214     tree_in -> SetBranchAddress("ClusPhiSize", &ClusPhiSize);
0215     tree_in -> SetBranchAddress("ClusZSize", &ClusZSize);
0216     
0217     tree_in -> SetBranchAddress("ClusLayer", &ClusLayer);
0218     tree_in -> SetBranchAddress("ClusLadderZId", &ClusLadderZId);
0219     tree_in -> SetBranchAddress("ClusLadderPhiId", &ClusLadderPhiId);
0220 
0221     tree_in -> SetBranchAddress("ClusEta_INTTz", &ClusEta_INTTz);
0222     tree_in -> SetBranchAddress("ClusEta_MBDz", &ClusEta_MBDz);
0223     tree_in -> SetBranchAddress("ClusPhi_AvgPV", &ClusPhi_AvgPV);
0224 
0225     // note : MC
0226     if (branch_map.find("ClusEta_TrueXYZ") != branch_map.end()) {tree_in -> SetBranchAddress("ClusEta_TrueXYZ", &ClusEta_TrueXYZ);}
0227     if (branch_map.find("ClusPhi_TrueXY") != branch_map.end()) {tree_in -> SetBranchAddress("ClusPhi_TrueXY", &ClusPhi_TrueXY);}
0228 
0229     // note : for data
0230     if (branch_map.find("MBDNSg2") != branch_map.end()) {tree_in -> SetBranchAddress("MBDNSg2", &MBDNSg2);}
0231     if (branch_map.find("MBDNSg2_vtxZ10cm") != branch_map.end()) {tree_in -> SetBranchAddress("MBDNSg2_vtxZ10cm", &MBDNSg2_vtxZ10cm);}
0232     if (branch_map.find("MBDNSg2_vtxZ30cm") != branch_map.end()) {tree_in -> SetBranchAddress("MBDNSg2_vtxZ30cm", &MBDNSg2_vtxZ30cm);}
0233 
0234     if (branch_map.find("InttBcoFullDiff_next") != branch_map.end()) {tree_in -> SetBranchAddress("InttBcoFullDiff_next", &InttBcoFullDiff_next); }
0235 
0236     // note : for MC
0237     if (branch_map.find("NTruthVtx") != branch_map.end()) {tree_in -> SetBranchAddress("NTruthVtx", &NTruthVtx);}
0238     if (branch_map.find("TruthPV_trig_x") != branch_map.end()) {tree_in -> SetBranchAddress("TruthPV_trig_x", &TruthPV_trig_x);}
0239     if (branch_map.find("TruthPV_trig_y") != branch_map.end()) {tree_in -> SetBranchAddress("TruthPV_trig_y", &TruthPV_trig_y);}
0240     if (branch_map.find("TruthPV_trig_z") != branch_map.end()) {tree_in -> SetBranchAddress("TruthPV_trig_z", &TruthPV_trig_z);}
0241 }
0242 
0243 void RestDist::PrepareOutputFile()
0244 {
0245     file_out = new TFile(Form("%s/%s", output_directory.c_str(), output_filename.c_str()), "RECREATE");
0246 }
0247 
0248 void RestDist::PrepareHist()
0249 {
0250     h1D_centrality_bin = new TH1D("h1D_centrality_bin","h1D_centrality_bin;Centrality [%];Entries",nCentrality_bin,&centrality_edges[0]); // note : the 0-5%
0251     
0252     h1D_map.clear();
0253     
0254     h1D_map.insert(
0255         std::make_pair(
0256             "h1D_PairDeltaEta",
0257             new TH1D("h1D_PairDeltaEta","h1D_PairDeltaEta;Pair #Delta#eta;Entries",nDeltaEtaBin,DeltaEtaEdge_min,DeltaEtaEdge_max)
0258         )
0259     );
0260 
0261     h1D_map.insert(
0262         std::make_pair(
0263             "h1D_PairDeltaPhi",
0264             new TH1D("h1D_PairDeltaPhi","h1D_PairDeltaPhi;Pair #Delta#phi;Entries",nDeltaPhiBin,DeltaPhiEdge_min,DeltaPhiEdge_max)
0265         )
0266     );
0267 
0268     h1D_map.insert(
0269         std::make_pair(
0270             "h1D_PairDeltaPhi_Narrow",
0271             new TH1D("h1D_PairDeltaPhi_Narrow","h1D_PairDeltaPhi_Narrow;Pair #Delta#phi;Entries",100,-0.05,0.05)
0272         )
0273     );
0274 
0275     h1D_map.insert(
0276         std::make_pair(
0277             "h1D_PairDeltaR",
0278             new TH1D("h1D_PairDeltaR","h1D_PairDeltaR;Pair #DeltaR;Entries",100, 0, 4)
0279         )
0280     );
0281 
0282     h1D_map.insert(
0283         std::make_pair(
0284             "h1D_PairDeltaR_Narrow",
0285             new TH1D("h1D_PairDeltaR_Narrow","h1D_PairDeltaR_Narrow;Pair #DeltaR;Entries",100, 0, 0.4)
0286         )
0287     );
0288 
0289     h1D_map.insert(
0290         std::make_pair(
0291             "h1D_PairDCA3D",
0292             new TH1D("h1D_PairDCA3D","h1D_PairDCA3D;Pair DCA 3D [cm];Entries",100, 0, 25)
0293         )
0294     );
0295 
0296     h1D_map.insert(std::make_pair("h1D_ClusPhi", new TH1D("h1D_ClusPhi", "h1D_ClusPhi;Cluster #phi (w.r.t beam spot);Entries (/0.05)",140,-3.5,3.5)));
0297     h1D_map.insert(std::make_pair("h1D_inner_ClusPhi", new TH1D("h1D_inner_ClusPhi", "h1D_inner_ClusPhi;Inner cluster #phi (w.r.t beam spot);Entries (/0.05)",140,-3.5,3.5)));
0298     h1D_map.insert(std::make_pair("h1D_outer_ClusPhi", new TH1D("h1D_outer_ClusPhi", "h1D_outer_ClusPhi;Outer cluster #phi (w.r.t beam spot);Entries (/0.05)",140,-3.5,3.5)));
0299 
0300     h1D_map.insert(std::make_pair("h1D_Cluster_phi_size", new TH1D("h1D_Cluster_phi_size", "h1D_Cluster_phi_size;Cluster #phi size;Entries (/1)",128,0,128)));
0301     h1D_map.insert(std::make_pair("h1D_Cluster_adc_size1", new TH1D("h1D_Cluster_adc_size1", "h1D_Cluster_adc_size1;Cluster Adc (PhiSize=1);Entries (/7)",25,0,250)));
0302     h1D_map.insert(std::make_pair("h1D_Cluster_adc", new TH1D("h1D_Cluster_adc", "h1D_Cluster_adc;Cluster Adc;Entries (/7)",100,0,700)));
0303     h1D_map.insert(std::make_pair("h1D_Cluster_Z", new TH1D("h1D_Cluster_Z", "h1D_Cluster_Z;Cluster Z [cm];Entries (/0.6)",100,-30,30)));
0304     h1D_map.insert(std::make_pair("h1D_Cluster_Z_typeA", new TH1D("h1D_Cluster_Z_typeA", "h1D_Cluster_Z_typeA;Cluster Z (Type A) [cm];Entries (/0.6)",100,-30,30)));
0305     h1D_map.insert(std::make_pair("h1D_Cluster_Z_typeB", new TH1D("h1D_Cluster_Z_typeB", "h1D_Cluster_Z_typeB;Cluster Z (Type B) [cm];Entries (/0.6)",100,-30,30)));
0306 
0307 
0308     h1D_map.insert(std::make_pair("h1D_ClusEta_INTTz", new TH1D("h1D_ClusEta_INTTz", "h1D_ClusEta_INTTz;Cluster #eta (INTTz);Entries (/0.1)",60,-3,3)));
0309     h1D_map.insert(std::make_pair("h1D_ClusEta_INTTz_highNClus", new TH1D("h1D_ClusEta_INTTz_highNClus", "h1D_ClusEta_INTTz_highNClus;Cluster #eta (INTTz);Entries (/0.1)",60,-3,3)));
0310     
0311     h1D_map.insert(std::make_pair("h1D_ClusEta_INTTz_typeA", new TH1D("h1D_ClusEta_INTTz_typeA", "h1D_ClusEta_INTTz_typeA;Cluster #eta (INTTz, typeA);Entries (/0.1)",60,-3,3)));
0312     h1D_map.insert(std::make_pair("h1D_ClusEta_INTTz_typeA_highNClus", new TH1D("h1D_ClusEta_INTTz_typeA_highNClus", "h1D_ClusEta_INTTz_typeA_highNClus;Cluster #eta (INTTz, typeA);Entries (/0.1)",60,-3,3)));
0313     
0314     h1D_map.insert(std::make_pair("h1D_ClusEta_INTTz_typeB", new TH1D("h1D_ClusEta_INTTz_typeB", "h1D_ClusEta_INTTz_typeB;Cluster #eta (INTTz, typeB);Entries (/0.1)",60,-3,3)));
0315     h1D_map.insert(std::make_pair("h1D_ClusEta_INTTz_typeB_highNClus", new TH1D("h1D_ClusEta_INTTz_typeB_highNClus", "h1D_ClusEta_INTTz_typeB_highNClus;Cluster #eta (INTTz, typeB);Entries (/0.1)",60,-3,3)));
0316     
0317     h1D_map.insert(std::make_pair("h1D_ClusEta_MBDz", new TH1D("h1D_ClusEta_MBDz", "h1D_ClusEta_MBDz;Cluster #eta (MBDz);Entries (/0.1)",60,-3,3)));
0318     h1D_map.insert(std::make_pair("h1D_ClusEta_MBDz_highNClus", new TH1D("h1D_ClusEta_MBDz_highNClus", "h1D_ClusEta_MBDz_highNClus;Cluster #eta (MBDz);Entries (/0.1)",60,-3,3)));
0319 
0320     h1D_map.insert(std::make_pair("h1D_ClusEta_MBDz_typeA", new TH1D("h1D_ClusEta_MBDz_typeA", "h1D_ClusEta_MBDz_typeA;Cluster #eta (MBDz, typeA);Entries (/0.1)",60,-3,3)));
0321     h1D_map.insert(std::make_pair("h1D_ClusEta_MBDz_typeA_highNClus", new TH1D("h1D_ClusEta_MBDz_typeA_highNClus", "h1D_ClusEta_MBDz_typeA_highNClus;Cluster #eta (MBDz, typeA);Entries (/0.1)",60,-3,3)));
0322     
0323     h1D_map.insert(std::make_pair("h1D_ClusEta_MBDz_typeB", new TH1D("h1D_ClusEta_MBDz_typeB", "h1D_ClusEta_MBDz_typeB;Cluster #eta (MBDz, typeB);Entries (/0.1)",60,-3,3)));
0324     h1D_map.insert(std::make_pair("h1D_ClusEta_MBDz_typeB_highNClus", new TH1D("h1D_ClusEta_MBDz_typeB_highNClus", "h1D_ClusEta_MBDz_typeB_highNClus;Cluster #eta (MBDz, typeB);Entries (/0.1)",60,-3,3)));
0325 
0326     if (runnumber == -1){
0327         h1D_map.insert(std::make_pair("h1D_ClusPhi_TrueXY", new TH1D("h1D_ClusPhi_TrueXY", "h1D_ClusPhi_TrueXY;Cluster #phi (w.r.t TrueXY);Entries (/0.05)",140,-3.5,3.5)));
0328 
0329         h1D_map.insert(std::make_pair("h1D_ClusEta_TrueXYZ", new TH1D("h1D_ClusEta_TrueXYZ", "h1D_ClusEta_TrueXYZ;Cluster #eta (TrueXYZ);Entries (/0.1)",60,-3,3)));
0330         h1D_map.insert(std::make_pair("h1D_ClusEta_TrueXYZ_highNClus", new TH1D("h1D_ClusEta_TrueXYZ_highNClus", "h1D_ClusEta_TrueXYZ_highNClus;Cluster #eta (TrueXYZ);Entries (/0.1)",60,-3,3)));
0331 
0332         h1D_map.insert(std::make_pair("h1D_ClusEta_TrueXYZ_typeA", new TH1D("h1D_ClusEta_TrueXYZ_typeA", "h1D_ClusEta_TrueXYZ_typeA;Cluster #eta (TrueXYZ, typeA);Entries (/0.1)",60,-3,3)));
0333         h1D_map.insert(std::make_pair("h1D_ClusEta_TrueXYZ_typeA_highNClus", new TH1D("h1D_ClusEta_TrueXYZ_typeA_highNClus", "h1D_ClusEta_TrueXYZ_typeA_highNClus;Cluster #eta (TrueXYZ, typeA);Entries (/0.1)",60,-3,3)));
0334         
0335         h1D_map.insert(std::make_pair("h1D_ClusEta_TrueXYZ_typeB", new TH1D("h1D_ClusEta_TrueXYZ_typeB", "h1D_ClusEta_TrueXYZ_typeB;Cluster #eta (TrueXYZ, typeB);Entries (/0.1)",60,-3,3)));
0336         h1D_map.insert(std::make_pair("h1D_ClusEta_TrueXYZ_typeB_highNClus", new TH1D("h1D_ClusEta_TrueXYZ_typeB_highNClus", "h1D_ClusEta_TrueXYZ_typeB_highNClus;Cluster #eta (TrueXYZ, typeB);Entries (/0.1)",60,-3,3)));
0337 
0338     }
0339 
0340     
0341     h1D_map.insert(std::make_pair("h1D_NClus", new TH1D("h1D_NClus", "h1D_NClus;NClus;Entries (/80)",100,0,10000)));
0342     h1D_map.insert(std::make_pair("h1D_NClus_inner", new TH1D("h1D_NClus_inner", "h1D_NClus_inner;NClus (inner);Entries (/45)",100,0,4500)));
0343     h1D_map.insert(std::make_pair("h1D_NClus_outer", new TH1D("h1D_NClus_outer", "h1D_NClus_outer;NClus (outer);Entries (/45)",100,0,4500)));
0344     
0345     h1D_map.insert(std::make_pair("h1D_inner_typeA_NClus", new TH1D("h1D_inner_typeA_NClus", "h1D_inner_typeA_NClus;NClus (inner, type A);Entries (/32)",100,0,3200)));
0346     h1D_map.insert(std::make_pair("h1D_outer_typeA_NClus", new TH1D("h1D_outer_typeA_NClus", "h1D_outer_typeA_NClus;NClus (outer, type A);Entries (/32)",100,0,3200)));
0347     h1D_map.insert(std::make_pair("h1D_inner_typeB_NClus", new TH1D("h1D_inner_typeB_NClus", "h1D_inner_typeB_NClus;NClus (inner, type B);Entries (/32)",100,0,3200)));
0348     h1D_map.insert(std::make_pair("h1D_outer_typeB_NClus", new TH1D("h1D_outer_typeB_NClus", "h1D_outer_typeB_NClus;NClus (outer, type B);Entries (/32)",100,0,3200)));
0349 
0350     h1D_map.insert(
0351         std::make_pair(
0352             "h1D_confirm_INTTvtxZ_Inclusive70",
0353             new TH1D("h1D_confirm_INTTvtxZ_Inclusive70","h1D_confirm_INTTvtxZ_Inclusive70;INTT vtx Z [cm];Entries (/0.6)",60,-60,60)
0354         )
0355     );
0356 
0357     h1D_map.insert(
0358         std::make_pair(
0359             "h1D_confirm_MBDvtxZ_Inclusive70",
0360             new TH1D("h1D_confirm_MBDvtxZ_Inclusive70","h1D_confirm_MBDvtxZ_Inclusive70;MBD vtx Z [cm];Entries (/0.6)",60,-60,60)
0361         )
0362     );
0363 
0364     for (auto &hist : h1D_map)
0365     {
0366         std::string YTitle = hist.second -> GetYaxis() -> GetTitle();
0367         std::string XTitle = hist.second -> GetXaxis() -> GetTitle();
0368 
0369         std::string YTitle_post;
0370 
0371         if (YTitle.find("Entries") != std::string::npos) // note : found the (Entries)
0372         {
0373             YTitle_post = Form("Entries  (/%.2f)", hist.second -> GetBinWidth(1));
0374             hist.second -> GetYaxis() -> SetTitle(YTitle_post.c_str());
0375         }
0376     }
0377 
0378     std::vector<double> used_edep_bins_vec = (runnumber != -1) ? data_edep_bins_vec : MC_edep_bins_vec;
0379     h1D_map.insert(std::make_pair("h1D_Cluster_adc_size1_real", new TH1D("h1D_Cluster_adc_size1_real", "h1D_Cluster_adc_size1_real;Cluster Adc (PhiSize=1);Entries",used_edep_bins_vec.size()-1,&used_edep_bins_vec[0])));
0380 
0381     // Division : ---h2D-----------------------------------------------------------------------------------------------
0382     h2D_map.clear();
0383 
0384     h2D_map.insert(
0385         std::make_pair(
0386             "h2D_confirm_InttNClus_MbdChargeSum",
0387             new TH2D("h2D_confirm_InttNClus_MbdChargeSum","h2D_confirm_InttNClus_MbdChargeSum;INTT NClus;MBD charge sum", 200, 0, 10000, 200, 0, 2500)
0388         )
0389     );
0390 
0391     h2D_map.insert(
0392         std::make_pair(
0393             "h2D_confirm_InttInnerOuterNClus",
0394             new TH2D("h2D_confirm_InttInnerOuterNClus","h2D_confirm_InttInnerOuterNClus;NClus (inner);NClus (outer)", 200, 0, 5500, 200, 0, 5500)
0395         )
0396     );
0397 
0398     h2D_map.insert(
0399         std::make_pair(
0400             "h2D_PairDCA3D_PairDeltaR",
0401             new TH2D("h2D_PairDCA3D_PairDeltaR","h2D_PairDCA3D_PairDeltaR;Pair DCA 3D [cm];Pair #DeltaR", 200, 0, 25, 200, 0, 4)
0402         )
0403     );
0404 
0405     h2D_map.insert(
0406         std::make_pair(
0407             "h2D_PairDCA3D_PairDeltaRNarrow",
0408             new TH2D("h2D_PairDCA3D_PairDeltaRNarrow","h2D_PairDCA3D_PairDeltaRNarrow;Pair DCA 3D [cm];Pair #DeltaR", 200, 0, 25, 200, 0, 0.4)
0409         )
0410     );
0411 
0412     h2D_map.insert(
0413         std::make_pair(
0414             "h2D_PairDCA3D_PairDeltaEta",
0415             new TH2D("h2D_PairDCA3D_PairDeltaEta","h2D_PairDCA3D_PairDeltaEta;Pair DCA 3D [cm];Pair #Delta#eta", 200, 0, 25, 200, -2, 2)
0416         )
0417     );
0418 
0419     h2D_map.insert(
0420         std::make_pair(
0421             "h2D_PairDCA3D_PairDeltaPhi",
0422             new TH2D("h2D_PairDCA3D_PairDeltaPhi","h2D_PairDCA3D_PairDeltaPhi;Pair DCA 3D [cm];Pair #Delta#phi", 200, 0, 25, nDeltaPhiBin, DeltaPhiEdge_min, DeltaPhiEdge_max)
0423         )
0424     );
0425 
0426     h2D_map.insert(
0427         std::make_pair(
0428             "h2D_TrackletEta_PairDeltaPhi",
0429             new TH2D("h2D_TrackletEta_PairDeltaPhi","h2D_TrackletEta_PairDeltaPhi;Tracklet #eta;Pair #Delta#phi", nEtaBin, EtaEdge_min, EtaEdge_max, nDeltaPhiBin, DeltaPhiEdge_min, DeltaPhiEdge_max)
0430         )
0431     );
0432 
0433     h2D_map.insert(
0434         std::make_pair(
0435             "h2D_PairDeltaPhi_DCAXY",
0436             new TH2D("h2D_PairDeltaPhi_DCAXY","h2D_PairDeltaPhi_DCAXY;Pair #Delta#phi;DCA 2D [cm]", nDeltaPhiBin, DeltaPhiEdge_min, DeltaPhiEdge_max, 100, -3, 3)
0437         )
0438     );
0439 
0440 
0441     h2D_map.insert(std::make_pair("h2D_ClusPhiSize_ClusAdc", new TH2D("h2D_ClusPhiSize_ClusAdc","h2D_ClusPhiSize_ClusAdc;Cluster #phi size;Cluster Adc",128,0,128,100,0,10000)));
0442     h2D_map.insert(std::make_pair("h2D_ClusXY", new TH2D("h2D_ClusXY","h2D_ClusXY;X [cm];Y [cm]",300,-15,15,300,-15,15)));
0443 
0444     // note : INTTz
0445     h2D_map.insert(std::make_pair("h2D_INTT_ClusEta_ClusAdc", new TH2D("h2D_INTT_ClusEta_ClusAdc","h2D_INTT_ClusEta_ClusAdc;Cluster #eta (INTTz);Cluster Adc",60,-3,3,50,0,500)));
0446     h2D_map.insert(std::make_pair("h2D_INTT_ClusEta_ClusAdc_highNClus", new TH2D("h2D_INTT_ClusEta_ClusAdc_highNClus","h2D_INTT_ClusEta_ClusAdc_highNClus;Cluster #eta (INTTz);Cluster Adc",60,-3,3,50,0,500)));
0447 
0448     // h2D_map.insert(std::make_pair("h2D_INTT_MBD_multiplicity", new TH2D("h2D_INTT_MBD_multiplicity","h2D_INTT_MBD_multiplicity;INTT NClus;MBD charge sum",200,0,10000,200,0,5500)));
0449     // h2D_map.insert(std::make_pair("h2D_inner_outer_INTT_multiplicity", new TH2D("h2D_inner_outer_INTT_multiplicity","h2D_inner_outer_INTT_multiplicity;NClus (inner);NClus (outer)",200,0,4000,200,0,4000)));
0450     h2D_map.insert(std::make_pair("h2D_ClusZ_LadderZID", new TH2D("h2D_ClusZ_LadderZID","h2D_ClusZ_LadderZID;Cluster Z [cm]; Clus Sensor ZId",100,-30,30,5,0,5)));
0451 
0452     h2D_map.insert(std::make_pair("h2D_sensor_HitMap", new TH2D("h2D_sensor_HitMap","h2D_sensor_HitMap;LayerID * 10 + sensor ZId;Ladder Phi Id",40,30,70,18,0,18)));
0453     h2D_map.insert(std::make_pair("h2D_sensor_AdcMap", new TH2D("h2D_sensor_AdcMap","h2D_sensor_AdcMap;LayerID * 10 + sensor ZId;Ladder Phi Id",40,30,70,18,0,18)));
0454     h2D_map.insert(std::make_pair("h2D_sensor_AvgAdcMap", new TH2D("h2D_sensor_AvgAdcMap","h2D_sensor_AvgAdcMap;LayerID * 10 + sensor ZId;Ladder Phi Id",40,30,70,18,0,18)));
0455 
0456     h2D_map.insert(std::make_pair("h2D_ClusXY_LargePhiSize", new TH2D("h2D_ClusXY_LargePhiSize","h2D_ClusXY_LargePhiSize;X [cm];Y [cm]",300,-15,15,300,-15,15)));    
0457     h2D_map.insert(std::make_pair("h2D_ClusXY_PhiSize43_46", new TH2D("h2D_ClusXY_PhiSize43_46","h2D_ClusXY_PhiSize43_46;X [cm];Y [cm]",300,-15,15,300,-15,15)));
0458     h2D_map.insert(std::make_pair("h2D_ClusXY_LargerPhiSize49", new TH2D("h2D_ClusXY_LargerPhiSize49","h2D_ClusXY_LargerPhiSize49;X [cm];Y [cm]",300,-15,15,300,-15,15)));    
0459 
0460     h2D_map.insert(std::make_pair("h2D_MBDcharge_south_north",   new TH2D("h2D_MBDcharge_south_north","h2D_MBDcharge_south_north;MBD charge sum (north);MBD charge sum (south)",100,0,1500,100,0,1500)));
0461     h2D_map.insert(std::make_pair("h2D_MBDvtxZ_MBDchargeSum",   new TH2D("h2D_MBDvtxZ_MBDchargeSum","h2D_MBDvtxZ_MBDchargeSum;MBD vtxZ [cm];MBD charge sum (all)",120,-60,60,100,0,3000)));
0462     h2D_map.insert(std::make_pair("h2D_MBDvtxZ_MBDchargeSouth", new TH2D("h2D_MBDvtxZ_MBDchargeSouth","h2D_MBDvtxZ_MBDchargeSouth;MBD vtxZ [cm];MBD charge sum (south)",120,-60,60,100,0,1500)));
0463     h2D_map.insert(std::make_pair("h2D_MBDvtxZ_MBDchargeNorth", new TH2D("h2D_MBDvtxZ_MBDchargeNorth","h2D_MBDvtxZ_MBDchargeNorth;MBD vtxZ [cm];MBD charge sum (north)",120,-60,60,100,0,1500)));
0464     h2D_map.insert(std::make_pair("h2D_MBD_vtxZ_ChargeAssy", new TH2D("h2D_MBD_vtxZ_ChargeAssy","h2D_MBD_vtxZ_ChargeAssy;MBD vtxZ [cm];MBD charge asymmetry",120,-60,60,100,-1.5,1.5)));
0465 
0466     h2D_map.insert(std::make_pair("h2D_NClus_ClusPhiSize", new TH2D("h2D_NClus_ClusPhiSize","h2D_NClus_ClusPhiSize;NClus;Cluster #phi size",100,0,10000,120,0,120)));
0467 }
0468 
0469 void RestDist::EvtCleanUp()
0470 {
0471     evt_sPH_inner_nocolumn_vec.clear();
0472     evt_sPH_outer_nocolumn_vec.clear();
0473 
0474     NClus_inner_typeA = 0;
0475     NClus_inner_typeB = 0;
0476     NClus_outer_typeA = 0;
0477     NClus_outer_typeB = 0;
0478 }
0479 
0480 void RestDist::PrepareClusterVec()
0481 {
0482     evt_sPH_inner_nocolumn_vec.clear();
0483     evt_sPH_outer_nocolumn_vec.clear();
0484 
0485     for (int clu_i = 0; clu_i < ClusX -> size(); clu_i++)
0486     {
0487         RestDist::clu_info this_clu;
0488 
0489         this_clu.adc = ClusAdc -> at(clu_i);
0490         this_clu.phi_size = ClusPhiSize -> at(clu_i);
0491         this_clu.sensorZID = ClusLadderZId -> at(clu_i);
0492         this_clu.ladderPhiID = ClusLadderPhiId -> at(clu_i);
0493         this_clu.layerID = ClusLayer -> at(clu_i);
0494 
0495         this_clu.index = clu_i;
0496 
0497         this_clu.x = ClusX -> at(clu_i);
0498         this_clu.y = ClusY -> at(clu_i);
0499         this_clu.z = ClusZ -> at(clu_i);
0500 
0501         this_clu.eta_INTTz = ClusEta_INTTz -> at(clu_i);
0502         this_clu.eta_MBDz  = ClusEta_MBDz  -> at(clu_i);
0503         this_clu.phi_avgXY = ClusPhi_AvgPV -> at(clu_i);
0504 
0505         if (runnumber == -1){
0506             this_clu.eta_TrueXYZ = ClusEta_TrueXYZ -> at(clu_i);
0507             this_clu.phi_TrueXY  = ClusPhi_TrueXY  -> at(clu_i);
0508         }
0509 
0510         if (isClusQA.first && this_clu.adc <= isClusQA.second.first) {continue;} // note : adc
0511         if (isClusQA.first && this_clu.phi_size > isClusQA.second.second) {continue;} // note : phi size
0512 
0513 
0514         if (this_clu.layerID == 3 || this_clu.layerID == 4){ // note : inner
0515             if (this_clu.sensorZID == typeA_sensorZID[0] || this_clu.sensorZID == typeA_sensorZID[1]){ // note : type A
0516                 NClus_inner_typeA++;
0517             }
0518             else { // note : type B
0519                 NClus_inner_typeB++;
0520             }
0521         }
0522         else { // note : outer 
0523             if (this_clu.sensorZID == typeA_sensorZID[0] || this_clu.sensorZID == typeA_sensorZID[1]){ // note : type A
0524                 NClus_outer_typeA++;
0525             }
0526             else { // note : type B
0527                 NClus_outer_typeB++;
0528             }
0529         }
0530 
0531         std::vector<RestDist::clu_info>* p_evt_sPH_nocolumn_vec =
0532         (this_clu.layerID == 3 || this_clu.layerID == 4) ? (&evt_sPH_inner_nocolumn_vec) : (&evt_sPH_outer_nocolumn_vec);
0533         p_evt_sPH_nocolumn_vec -> push_back(this_clu);
0534     }
0535 }
0536 
0537 void RestDist::PrepareEvent()
0538 {
0539     for (int i = 0; i < run_nEvents; i++)
0540     {
0541         tree_in -> GetEntry(i);
0542 
0543         EvtCleanUp();
0544 
0545         if (i % 10 == 0) {std::cout << "Processing event " << i << std::endl;}
0546 
0547         // =======================================================================================================================================================
0548         // note : hard cut
0549 
0550         // note : for data
0551         if (runnumber != -1 && ApplyEvtBcoFullDiffCut.first && InttBcoFullDiff_next <= ApplyEvtBcoFullDiffCut.second) {continue;}
0552         if (runnumber != -1 && MBDNSg2 != 1) {continue;} // todo: assume MC no trigger
0553 
0554         // note : for MC
0555         // if (runnumber == -1 && NTruthVtx != 1) {continue;}
0556 
0557         // note : both data and MC
0558         if (is_min_bias != 1) {continue;}
0559         if (MBD_z_vtx != MBD_z_vtx) {continue;}
0560         if (MBD_centrality != MBD_centrality) {continue;}
0561         if (MBD_centrality < 0 || MBD_centrality > 100) {continue;}
0562         if (INTTvtxZ != INTTvtxZ) {continue;}
0563         if (MBD_z_vtx < cut_GlobalMBDvtxZ.first || MBD_z_vtx > cut_GlobalMBDvtxZ.second) {continue;} // todo: the hard cut 60 cm 
0564 
0565         // =======================================================================================================================================================
0566         // note : optional cut
0567         if (Apply_cut && (MBD_z_vtx - INTTvtxZ < cut_vtxZDiff.first || MBD_z_vtx - INTTvtxZ > cut_vtxZDiff.second) ) {continue;}
0568         if (Apply_cut && (TrapezoidalFitWidth < cut_TrapezoidalFitWidth.first || TrapezoidalFitWidth > cut_TrapezoidalFitWidth.second)){continue;}
0569         if (Apply_cut && (TrapezoidalFWHM < cut_TrapezoidalFWHM.first || TrapezoidalFWHM > cut_TrapezoidalFWHM.second)){continue;}
0570         if (Apply_cut && (INTTvtxZError < cut_INTTvtxZError.first || INTTvtxZError > cut_INTTvtxZError.second)){continue;}
0571 
0572         // =======================================================================================================================================================
0573         if (ApplyVtxZReWeighting && runnumber != -1){
0574             std::cout<<"Should not have the vtxZ weighting from the data"<<std::endl;
0575             exit(1);
0576         }
0577 
0578         double INTTvtxZWeighting;
0579         if (ApplyVtxZReWeighting && h1D_INTT_vtxZ_reweighting != nullptr){
0580             INTTvtxZWeighting = h1D_INTT_vtxZ_reweighting -> GetBinContent(h1D_INTT_vtxZ_reweighting -> FindBin(INTTvtxZ));
0581         }
0582         else if (ApplyVtxZReWeighting && h1D_INTT_vtxZ_reweighting == nullptr){
0583             std::cout << "ApplyVtxZReWeighting is true, but h1D_INTT_vtxZ_reweighting is nullptr" << std::endl;
0584             exit(1);
0585         }
0586         else {
0587             INTTvtxZWeighting = 1.0;
0588         }
0589         
0590         int Mbin = h1D_centrality_bin -> Fill(MBD_centrality, INTTvtxZWeighting);
0591         Mbin = (Mbin == -1) ? -1 : Mbin - 1;
0592         if (Mbin == -1) {
0593             std::cout << "Mbin == -1, MBD_centrality = " << MBD_centrality << std::endl;
0594             continue;
0595         }
0596 
0597         // =======================================================================================================================================================
0598         // todo: Only use the INTT for the analysis vtxZ range
0599         if (RequireVtxZRange.first && (INTTvtxZ < RequireVtxZRange.second.first || INTTvtxZ > RequireVtxZRange.second.second)) {continue;}
0600         // if (RequireVtxZRange.first && (MBD_z_vtx < RequireVtxZRange.second.first || MBD_z_vtx > RequireVtxZRange.second.second)) {continue;}
0601          
0602         PrepareClusterVec();
0603 
0604         if(isRotated) {GetRotatedClusterVec(evt_sPH_inner_nocolumn_vec);}
0605 
0606         if (Mbin <= Semi_inclusive_bin){
0607             h1D_map["h1D_confirm_INTTvtxZ_Inclusive70"] -> Fill(INTTvtxZ, INTTvtxZWeighting);
0608             h1D_map["h1D_confirm_MBDvtxZ_Inclusive70"] -> Fill(MBD_z_vtx, INTTvtxZWeighting);
0609         }
0610 
0611         int total_NClus = NClus_inner_typeA + NClus_inner_typeB + NClus_outer_typeA + NClus_outer_typeB;
0612 
0613         h1D_map["h1D_NClus"] -> Fill(total_NClus, INTTvtxZWeighting);
0614         h1D_map["h1D_NClus_inner"] -> Fill(NClus_inner_typeA + NClus_inner_typeB, INTTvtxZWeighting);
0615         h1D_map["h1D_NClus_outer"] -> Fill(NClus_outer_typeA + NClus_outer_typeB, INTTvtxZWeighting);
0616         h1D_map["h1D_inner_typeA_NClus"] -> Fill(NClus_inner_typeA, INTTvtxZWeighting);
0617         h1D_map["h1D_inner_typeB_NClus"] -> Fill(NClus_inner_typeB, INTTvtxZWeighting);
0618         h1D_map["h1D_outer_typeA_NClus"] -> Fill(NClus_outer_typeA, INTTvtxZWeighting);
0619         h1D_map["h1D_outer_typeB_NClus"] -> Fill(NClus_outer_typeB, INTTvtxZWeighting);
0620 
0621         // note : h2D
0622         h2D_map["h2D_confirm_InttNClus_MbdChargeSum"] -> Fill(total_NClus, MBD_charge_sum, INTTvtxZWeighting);
0623         h2D_map["h2D_confirm_InttInnerOuterNClus"] -> Fill(NClus_inner_typeA + NClus_inner_typeB, NClus_outer_typeA + NClus_outer_typeB, INTTvtxZWeighting);
0624         h2D_map["h2D_MBD_vtxZ_ChargeAssy"] -> Fill(MBD_z_vtx, MBD_charge_asymm, INTTvtxZWeighting);
0625         h2D_map["h2D_MBDcharge_south_north"] -> Fill(MBD_south_charge_sum, MBD_north_charge_sum, INTTvtxZWeighting);
0626         h2D_map["h2D_MBDvtxZ_MBDchargeSum"] -> Fill(MBD_z_vtx, MBD_charge_sum, INTTvtxZWeighting);
0627         h2D_map["h2D_MBDvtxZ_MBDchargeSouth"] -> Fill(MBD_z_vtx, MBD_south_charge_sum, INTTvtxZWeighting);
0628         h2D_map["h2D_MBDvtxZ_MBDchargeNorth"] -> Fill(MBD_z_vtx, MBD_north_charge_sum, INTTvtxZWeighting);
0629 
0630         for (int layer_index = 0; layer_index < 2; layer_index++){ // note : inner and outer barrels
0631             std::vector<RestDist::clu_info>* p_evt_sPH_nocolumn_vec = (layer_index == 0) ? (&evt_sPH_inner_nocolumn_vec) : (&evt_sPH_outer_nocolumn_vec);
0632 
0633             for (int clu_i = 0; clu_i < p_evt_sPH_nocolumn_vec->size(); clu_i++){
0634                 
0635                 RestDist::clu_info this_clu = p_evt_sPH_nocolumn_vec -> at(clu_i);
0636 
0637                 h1D_map["h1D_Cluster_Z"] -> Fill(this_clu.z, INTTvtxZWeighting);
0638                 h1D_map["h1D_ClusEta_MBDz"] -> Fill(this_clu.eta_MBDz, INTTvtxZWeighting);
0639                 h1D_map["h1D_ClusPhi"] -> Fill(this_clu.phi_avgXY, INTTvtxZWeighting);
0640                 h1D_map["h1D_ClusEta_INTTz"] -> Fill(this_clu.eta_INTTz, INTTvtxZWeighting);
0641                 h1D_map["h1D_Cluster_phi_size"] -> Fill(this_clu.phi_size, INTTvtxZWeighting);
0642                 h1D_map["h1D_Cluster_adc"] -> Fill(this_clu.adc, INTTvtxZWeighting);
0643 
0644                 if (layer_index == 0){
0645                     h1D_map["h1D_inner_ClusPhi"] -> Fill(this_clu.phi_avgXY, INTTvtxZWeighting);
0646                 }
0647                 else if (layer_index == 1){
0648                     h1D_map["h1D_outer_ClusPhi"] -> Fill(this_clu.phi_avgXY, INTTvtxZWeighting);
0649                 }
0650                 
0651                 if (this_clu.phi_size == 1) {h1D_map["h1D_Cluster_adc_size1"] -> Fill(this_clu.adc, INTTvtxZWeighting);}
0652                 if (this_clu.phi_size == 1) {h1D_map["h1D_Cluster_adc_size1_real"] -> Fill(this_clu.adc, INTTvtxZWeighting);}
0653 
0654                 if (runnumber == -1) {h1D_map["h1D_ClusPhi_TrueXY"] -> Fill(this_clu.phi_TrueXY), INTTvtxZWeighting;}
0655                 if (runnumber == -1) {h1D_map["h1D_ClusEta_TrueXYZ"] -> Fill(this_clu.eta_TrueXYZ), INTTvtxZWeighting;}
0656 
0657                 
0658 
0659                 if (total_NClus > HighNClus){
0660                     h1D_map["h1D_ClusEta_MBDz_highNClus"] -> Fill(this_clu.eta_MBDz, INTTvtxZWeighting);
0661                     h1D_map["h1D_ClusEta_INTTz_highNClus"] -> Fill(this_clu.eta_INTTz, INTTvtxZWeighting);
0662 
0663                     if (runnumber == -1) { h1D_map["h1D_ClusEta_TrueXYZ_highNClus"] -> Fill(this_clu.eta_TrueXYZ), INTTvtxZWeighting;}
0664                 }
0665 
0666                 // note : h2D
0667                 h2D_map["h2D_ClusPhiSize_ClusAdc"] -> Fill(this_clu.phi_size, this_clu.adc, INTTvtxZWeighting);
0668                 h2D_map["h2D_ClusXY"] -> Fill(this_clu.x, this_clu.y, INTTvtxZWeighting);
0669                 h2D_map["h2D_INTT_ClusEta_ClusAdc"] -> Fill(this_clu.eta_INTTz, this_clu.adc, INTTvtxZWeighting);
0670                 h2D_map["h2D_ClusZ_LadderZID"] -> Fill(this_clu.z, this_clu.sensorZID, INTTvtxZWeighting);
0671                 h2D_map["h2D_sensor_HitMap"] -> Fill(this_clu.layerID * 10 + this_clu.sensorZID, this_clu.ladderPhiID, INTTvtxZWeighting);
0672                 h2D_map["h2D_sensor_AdcMap"] -> Fill(this_clu.layerID * 10 + this_clu.sensorZID, this_clu.ladderPhiID, this_clu.adc * INTTvtxZWeighting);
0673 
0674                 if (this_clu.phi_size > 30) {h2D_map["h2D_ClusXY_LargePhiSize"]->Fill(this_clu.x, this_clu.y, INTTvtxZWeighting);}
0675                 if (this_clu.phi_size == 43 || this_clu.phi_size == 46){h2D_map["h2D_ClusXY_PhiSize43_46"] -> Fill(this_clu.x, this_clu.y, INTTvtxZWeighting);}
0676                 if (this_clu.phi_size >= 49){h2D_map["h2D_ClusXY_LargerPhiSize49"]->Fill(this_clu.x, this_clu.y, INTTvtxZWeighting);}
0677                 h2D_map["h2D_NClus_ClusPhiSize"] -> Fill(total_NClus, this_clu.phi_size, INTTvtxZWeighting);
0678 
0679                 if (total_NClus > HighNClus){
0680                     h2D_map["h2D_INTT_ClusEta_ClusAdc_highNClus"] -> Fill(this_clu.eta_INTTz, this_clu.adc, INTTvtxZWeighting);
0681                 }
0682 
0683                 // Division : ---type--------------------------------------------------------------------------------------------------------------------------
0684                 if (this_clu.sensorZID == typeA_sensorZID[0] || this_clu.sensorZID == typeA_sensorZID[1]){ // note : type A
0685                     h1D_map["h1D_Cluster_Z_typeA"] -> Fill(this_clu.z, INTTvtxZWeighting);
0686                     h1D_map["h1D_ClusEta_MBDz_typeA"] -> Fill(this_clu.eta_MBDz, INTTvtxZWeighting);
0687                     h1D_map["h1D_ClusEta_INTTz_typeA"] -> Fill(this_clu.eta_INTTz, INTTvtxZWeighting);
0688 
0689                     if (runnumber == -1) {h1D_map["h1D_ClusEta_TrueXYZ_typeA"] -> Fill(this_clu.eta_TrueXYZ), INTTvtxZWeighting;}
0690 
0691                     if (total_NClus > HighNClus){
0692                         h1D_map["h1D_ClusEta_MBDz_typeA_highNClus"] -> Fill(this_clu.eta_MBDz, INTTvtxZWeighting);
0693                         h1D_map["h1D_ClusEta_INTTz_typeA_highNClus"] -> Fill(this_clu.eta_INTTz, INTTvtxZWeighting);
0694 
0695                         if (runnumber == -1) {h1D_map["h1D_ClusEta_TrueXYZ_typeA_highNClus"] -> Fill(this_clu.eta_TrueXYZ, INTTvtxZWeighting);}
0696                     }
0697                 }
0698                 else { // note : type B
0699                     h1D_map["h1D_Cluster_Z_typeB"] -> Fill(this_clu.z, INTTvtxZWeighting);
0700                     h1D_map["h1D_ClusEta_MBDz_typeB"] -> Fill(this_clu.eta_MBDz, INTTvtxZWeighting);
0701                     h1D_map["h1D_ClusEta_INTTz_typeB"] -> Fill(this_clu.eta_INTTz, INTTvtxZWeighting);
0702 
0703                     if (runnumber == -1) {h1D_map["h1D_ClusEta_TrueXYZ_typeB"] -> Fill(this_clu.eta_TrueXYZ), INTTvtxZWeighting;}
0704 
0705                     if (total_NClus > HighNClus){
0706                         h1D_map["h1D_ClusEta_MBDz_typeB_highNClus"] -> Fill(this_clu.eta_MBDz, INTTvtxZWeighting);
0707                         h1D_map["h1D_ClusEta_INTTz_typeB_highNClus"] -> Fill(this_clu.eta_INTTz, INTTvtxZWeighting);
0708 
0709                         if (runnumber == -1) {h1D_map["h1D_ClusEta_TrueXYZ_typeB_highNClus"] -> Fill(this_clu.eta_TrueXYZ, INTTvtxZWeighting);}
0710                     }
0711                 }
0712 
0713             }
0714         }
0715 
0716 
0717         // Division : --- For pairs ---------------------------------------------------------------------------------------------------------------------
0718         inner_clu_phi_map.clear();
0719         outer_clu_phi_map.clear();
0720         inner_clu_phi_map = std::vector<std::vector<std::pair<bool,RestDist::clu_info>>>(360);
0721         outer_clu_phi_map = std::vector<std::vector<std::pair<bool,RestDist::clu_info>>>(360);
0722 
0723         for (int inner_i = 0; inner_i < int(evt_sPH_inner_nocolumn_vec.size()); inner_i++) {
0724             RestDist::clu_info this_clu = evt_sPH_inner_nocolumn_vec[inner_i];
0725             double Clus_InnerPhi_Offset = (this_clu.phi_avgXY < 0) ? this_clu.phi_avgXY * (180./TMath::Pi()) + 360 : this_clu.phi_avgXY * (180./TMath::Pi());
0726             inner_clu_phi_map[ int(Clus_InnerPhi_Offset) ].push_back({false,this_clu});
0727         }
0728         for (int outer_i = 0; outer_i < int(evt_sPH_outer_nocolumn_vec.size()); outer_i++) {
0729             RestDist::clu_info this_clu = evt_sPH_outer_nocolumn_vec[outer_i];
0730             double Clus_OuterPhi_Offset = (this_clu.phi_avgXY < 0) ? this_clu.phi_avgXY * (180./TMath::Pi()) + 360 : this_clu.phi_avgXY * (180./TMath::Pi());
0731             outer_clu_phi_map[ int(Clus_OuterPhi_Offset) ].push_back({false,this_clu});            
0732         }
0733 
0734         // temp_all_pair_vec.clear();
0735         // temp_all_pair_deltaR_vec.clear();
0736 
0737         for (int inner_phi_i = 0; inner_phi_i < 360; inner_phi_i++) // note : each phi cell (1 degree)
0738         {
0739             // note : N cluster in this phi cell
0740             for (int inner_phi_clu_i = 0; inner_phi_clu_i < int(inner_clu_phi_map[inner_phi_i].size()); inner_phi_clu_i++)
0741             {
0742                 RestDist::clu_info inner_clu = inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second;
0743 
0744                 // todo: change the outer phi scan range
0745                 // note : the outer phi index, -1, 0, 1
0746                 // note : the outer phi index, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 for the scan test
0747                 for (int scan_i = -5; scan_i < 6; scan_i++)
0748                 {
0749                     int true_scan_i = ((inner_phi_i + scan_i) < 0) ? 360 + (inner_phi_i + scan_i) : ((inner_phi_i + scan_i) > 359) ? (inner_phi_i + scan_i)-360 : inner_phi_i + scan_i;
0750 
0751                     // note : N clusters in that outer phi cell
0752                     for (int outer_phi_clu_i = 0; outer_phi_clu_i < int(outer_clu_phi_map[true_scan_i].size()); outer_phi_clu_i++)
0753                     {                        
0754                         RestDist::clu_info outer_clu = outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second;
0755 
0756                         double PairEta = (inner_clu.eta_INTTz + outer_clu.eta_INTTz)/2.;
0757                         double deltaEta = inner_clu.eta_INTTz - outer_clu.eta_INTTz;
0758                         double deltaPhi = inner_clu.phi_avgXY - outer_clu.phi_avgXY;
0759                         double deltaR = sqrt(pow(deltaPhi,2)+pow(deltaEta,2));
0760                         double DCAXY_sign = calculateAngleBetweenVectors(
0761                             outer_clu.x, outer_clu.y,
0762                             inner_clu.x, inner_clu.y,
0763                             vertexXYIncm.first, vertexXYIncm.second
0764                         );
0765                         double PairRadius = shortestDistanceLineToPoint(
0766                             {inner_clu.x, inner_clu.y, inner_clu.z},
0767                             {outer_clu.x, outer_clu.y, outer_clu.z},
0768                             {vertexXYIncm.first, vertexXYIncm.second, INTTvtxZ}
0769                         );
0770 
0771                         h1D_map["h1D_PairDeltaEta"] -> Fill(deltaEta, INTTvtxZWeighting);
0772                         h1D_map["h1D_PairDeltaPhi"] -> Fill(deltaPhi, INTTvtxZWeighting);
0773                         h1D_map["h1D_PairDeltaPhi_Narrow"] -> Fill(deltaPhi, INTTvtxZWeighting);
0774 
0775                         h1D_map["h1D_PairDeltaR"] -> Fill(deltaR, INTTvtxZWeighting);
0776                         h1D_map["h1D_PairDeltaR_Narrow"] -> Fill(deltaR, INTTvtxZWeighting);
0777 
0778                         h1D_map["h1D_PairDCA3D"] -> Fill(PairRadius, INTTvtxZWeighting);
0779 
0780                         // todo: require the pair can be linked to the vertex
0781                         if (fabs(deltaEta) <= 0.25){h2D_map["h2D_TrackletEta_PairDeltaPhi"] -> Fill(PairEta,deltaPhi,INTTvtxZWeighting);}
0782                         h2D_map["h2D_PairDeltaPhi_DCAXY"] -> Fill(deltaPhi, DCAXY_sign,INTTvtxZWeighting);
0783 
0784                         h2D_map["h2D_PairDCA3D_PairDeltaR"] -> Fill(PairRadius,deltaR,INTTvtxZWeighting);
0785                         h2D_map["h2D_PairDCA3D_PairDeltaRNarrow"] -> Fill(PairRadius,deltaR,INTTvtxZWeighting);
0786                         h2D_map["h2D_PairDCA3D_PairDeltaEta"] -> Fill(PairRadius,deltaEta,INTTvtxZWeighting);
0787                         h2D_map["h2D_PairDCA3D_PairDeltaPhi"] -> Fill(PairRadius,deltaPhi,INTTvtxZWeighting);
0788 
0789                         if (i%50 == 0 && fabs(deltaPhi) <= 0.021 && fabs(deltaEta) <= 0.25){
0790                             std::cout<<"deltaPhi = "<<deltaPhi<<", deltaEta = "<<deltaEta<<std::endl;
0791                             std::cout<<"inner_clu : {"<<inner_clu.x<<", "<<inner_clu.y<<", "<<inner_clu.z<<"}"<<std::endl;
0792                             std::cout<<"outer_clu : {"<<outer_clu.x<<", "<<outer_clu.y<<", "<<outer_clu.z<<"}"<<std::endl;
0793                             std::cout<<"vertex : {"<<vertexXYIncm.first<<", "<<vertexXYIncm.second<<", "<<INTTvtxZ<<"}"<<std::endl;
0794                             std::cout<<"PairRadius = "<<PairRadius<<std::endl;
0795                             std::cout<<std::endl;
0796                         } 
0797 
0798 
0799                         // pair_str this_pair;
0800                         // this_pair.delta_phi = deltaPhi;
0801                         // this_pair.delta_eta = deltaEta;
0802                         // this_pair.inner_index = inner_clu.index;                        
0803                         // this_pair.outer_index = outer_clu.index;
0804                         // temp_all_pair_vec.push_back(this_pair);
0805                         // temp_all_pair_deltaR_vec.push_back(deltaR);
0806 
0807                     }
0808                 }
0809             }
0810         }         
0811 
0812         // long long vec_size = temp_all_pair_deltaR_vec.size();
0813         // long long ind[temp_all_pair_deltaR_vec.size()];
0814         // TMath::Sort(vec_size, &temp_all_pair_deltaR_vec[0], ind, false);
0815 
0816         // BestPair_h2D_map.insert(
0817         //     std::make_pair(
0818         //         i,
0819         //         new TH2D(Form("h2D_BestPair_%d",i),Form("h2D_BestPair_%d;Pair Delta Eta;Pair Delta Phi",i),400,-0.5,0.5,400,-0.5,0.5)
0820         //     )
0821         // );
0822 
0823         // BestPair_gr_map.insert(
0824         //     std::make_pair(
0825         //         i,
0826         //         new TGraph()
0827         //     )
0828         // ); BestPair_gr_map[i] -> SetTitle(Form("h2D_gr_%d;pair_index;#Delta#phi [rad]",i));
0829 
0830         // std::vector<int> Used_Clus_index_vec; Used_Clus_index_vec.clear();
0831         // for (int pair_i = 0; pair_i < vec_size; pair_i++){
0832 
0833         //     pair_str pair = temp_all_pair_vec[ind[pair_i]];
0834         //     if (std::find(Used_Clus_index_vec.begin(), Used_Clus_index_vec.end(), pair.inner_index) != Used_Clus_index_vec.end()) {continue;}
0835         //     if (std::find(Used_Clus_index_vec.begin(), Used_Clus_index_vec.end(), pair.outer_index) != Used_Clus_index_vec.end()) {continue;}
0836 
0837         //     if (temp_all_pair_deltaR_vec[ind[pair_i]] > 0.5) {break;}
0838 
0839         //     BestPair_h2D_map[i] -> SetBinContent(
0840         //         BestPair_h2D_map[i] -> FindBin(pair.delta_eta, pair.delta_phi),
0841         //         pair_i
0842         //     );
0843 
0844         //     BestPair_gr_map[i] -> SetPoint(BestPair_gr_map[i] -> GetN(), pair_i, pair.delta_phi);
0845 
0846         //     Used_Clus_index_vec.push_back(pair.inner_index);
0847         //     Used_Clus_index_vec.push_back(pair.outer_index);
0848         // }
0849 
0850         // BestPair_gr_map[i]->GetYaxis()->SetRangeUser(-0.5,0.5);
0851 
0852 
0853         // for (int inner_i = 0; inner_i < evt_sPH_inner_nocolumn_vec.size(); inner_i++){
0854         //     for (int outer_i = 0; outer_i < evt_sPH_outer_nocolumn_vec.size(); outer_i++){
0855 
0856         //         RestDist::clu_info inner_clu = evt_sPH_inner_nocolumn_vec.at(inner_i);
0857         //         RestDist::clu_info outer_clu = evt_sPH_outer_nocolumn_vec.at(outer_i);
0858 
0859         //         double deltaEta = inner_clu.eta_INTTz - outer_clu.eta_INTTz;
0860         //         double deltaPhi = inner_clu.phi_avgXY - outer_clu.phi_avgXY;
0861         //         double deltaR = sqrt(pow(deltaPhi,2)+pow(deltaEta,2));
0862 
0863         //         h1D_map["h1D_PairDeltaEta"] -> Fill(deltaEta, INTTvtxZWeighting);
0864         //         h1D_map["h1D_PairDeltaPhi"] -> Fill(deltaPhi, INTTvtxZWeighting);
0865         //         h1D_map["h1D_PairDeltaPhi_Narrow"] -> Fill(deltaPhi, INTTvtxZWeighting);
0866 
0867         //         h1D_map["h1D_PairDeltaR"] -> Fill(deltaR, INTTvtxZWeighting);
0868 
0869         //         double PairRadius = shortestDistanceLineToPoint(
0870         //             {inner_clu.x, inner_clu.y, inner_clu.z},
0871         //             {outer_clu.x, outer_clu.y, outer_clu.z},
0872         //             {vertexXYIncm.first, vertexXYIncm.second, INTTvtxZ}
0873         //         );
0874         //         h1D_map["h1D_PairDCA3D"] -> Fill(PairRadius, INTTvtxZWeighting);
0875 
0876         //         h2D_map["h2D_PairDCA3D_PairDeltaR"] -> Fill(PairRadius,deltaR,INTTvtxZWeighting);
0877         //         h2D_map["h2D_PairDCA3D_PairDeltaEta"] -> Fill(PairRadius,deltaEta,INTTvtxZWeighting);
0878         //         h2D_map["h2D_PairDCA3D_PairDeltaPhi"] -> Fill(PairRadius,deltaPhi,INTTvtxZWeighting);
0879 
0880         //         if (i%50 == 0 && fabs(deltaPhi) <= 0.021 && fabs(deltaEta) <= 0.25){
0881         //             std::cout<<"deltaPhi = "<<deltaPhi<<", deltaEta = "<<deltaEta<<std::endl;
0882         //             std::cout<<"inner_clu : {"<<inner_clu.x<<", "<<inner_clu.y<<", "<<inner_clu.z<<"}"<<std::endl;
0883         //             std::cout<<"outer_clu : {"<<outer_clu.x<<", "<<outer_clu.y<<", "<<outer_clu.z<<"}"<<std::endl;
0884         //             std::cout<<"vertex : {"<<vertexXYIncm.first<<", "<<vertexXYIncm.second<<", "<<INTTvtxZ<<"}"<<std::endl;
0885         //             std::cout<<"PairRadius = "<<PairRadius<<std::endl;
0886         //             std::cout<<std::endl;
0887         //         } 
0888         //     }
0889         // }
0890         
0891     }// note : end of event loop
0892 
0893 }
0894 
0895 void RestDist::EndRun()
0896 {
0897     file_out -> cd();
0898     h1D_centrality_bin -> Write();
0899 
0900     // for (auto &pair : BestPair_h2D_map)
0901     // {
0902     //     pair.second -> Write();
0903     //     BestPair_gr_map[pair.first] -> Write(Form("h2D_gr_%d",pair.first));
0904     // }
0905 
0906     for (auto &hist : h1D_map)
0907     {
0908         hist.second -> Write();
0909     }
0910 
0911     for (auto &hist : h2D_map)
0912     {
0913         hist.second -> Write();
0914     }
0915 
0916     h2D_map["h2D_sensor_AdcMap"] -> Sumw2(true);
0917     h2D_map["h2D_sensor_HitMap"] -> Sumw2(true);
0918     h2D_map["h2D_sensor_AvgAdcMap"] -> Divide(h2D_map["h2D_sensor_AdcMap"], h2D_map["h2D_sensor_HitMap"], 1, 1);
0919     h2D_map["h2D_sensor_AvgAdcMap"] -> Write();
0920 
0921     file_out -> Close();
0922 }
0923 
0924 double RestDist::shortestDistanceLineToPoint(const std::vector<double> point1, const std::vector<double> point2, const std::vector<double> target) {
0925     // note : Ensure the vectors are 3-dimensional
0926     if (point1.size() != 3 || point2.size() != 3 || target.size() != 3) {
0927         throw std::invalid_argument("All input vectors must have exactly 3 elements.");
0928         exit(1);
0929     }
0930 
0931     // note : Compute direction vector of the line (point2 - point1)
0932     // note : from point1 to point2
0933     std::vector<double> lineDir = {
0934         point2[0] - point1[0],
0935         point2[1] - point1[1],
0936         point2[2] - point1[2]
0937     };
0938 
0939     // note : Compute vector from point1 to the target point
0940     // note : from point1 to target
0941     std::vector<double> pointToTarget = {
0942         target[0] - point1[0],
0943         target[1] - point1[1],
0944         target[2] - point1[2]
0945     };
0946 
0947     // note : Compute the cross product of lineDir and pointToTarget
0948     std::vector<double> crossProduct = {
0949         lineDir[1] * pointToTarget[2] - lineDir[2] * pointToTarget[1],
0950         lineDir[2] * pointToTarget[0] - lineDir[0] * pointToTarget[2],
0951         lineDir[0] * pointToTarget[1] - lineDir[1] * pointToTarget[0]
0952     };
0953 
0954     // note : Compute the magnitude of the cross product
0955     double crossProductMagnitude = std::sqrt( pow(crossProduct[0],2) + pow(crossProduct[1],2) + pow(crossProduct[2],2) );
0956 
0957     // note : Compute the magnitude of the line direction vector
0958     double lineDirMagnitude = std::sqrt( pow(lineDir[0],2) + pow(lineDir[1],2) + pow(lineDir[2],2) );
0959 
0960     // note : Compute and return the shortest distance
0961     return crossProductMagnitude / lineDirMagnitude;
0962 }
0963 
0964 void RestDist::GetRotatedClusterVec(std::vector<RestDist::clu_info> &input_cluster_vec)
0965 {
0966     for (RestDist::clu_info &this_clu : input_cluster_vec)
0967     {
0968         std::pair<double,double> rotated_xy = rotatePoint(this_clu.x, this_clu.y);
0969      
0970         this_clu.x = rotated_xy.first;
0971         this_clu.y = rotated_xy.second;
0972 
0973         this_clu.eta_INTTz = RestDist::get_clu_eta({vertexXYIncm.first, vertexXYIncm.second, INTTvtxZ}, {this_clu.x, this_clu.y, this_clu.z});
0974         this_clu.eta_MBDz = RestDist::get_clu_eta({vertexXYIncm.first, vertexXYIncm.second, MBD_z_vtx}, {this_clu.x, this_clu.y, this_clu.z});
0975         this_clu.phi_avgXY = atan2(this_clu.y - vertexXYIncm.second, this_clu.x - vertexXYIncm.first);
0976 
0977         if (runnumber == -1){
0978             this_clu.eta_TrueXYZ = RestDist::get_clu_eta({TruthPV_trig_x, TruthPV_trig_y, TruthPV_trig_z}, {this_clu.x, this_clu.y, this_clu.z});
0979             this_clu.phi_TrueXY = atan2(this_clu.y - TruthPV_trig_y, this_clu.x - TruthPV_trig_x);
0980         }
0981 
0982     }
0983 }
0984 
0985 std::pair<double,double> RestDist::rotatePoint(double x, double y)
0986 {
0987     // Convert the rotation angle from degrees to radians
0988     double rotation = rotate_phi_angle; // todo rotation is here
0989     double angleRad = rotation * M_PI / 180.0;
0990 
0991     // Perform the rotation
0992     double xOut = x * cos(angleRad) - y * sin(angleRad);
0993     double yOut = x * sin(angleRad) + y * cos(angleRad);
0994 
0995     xOut = (fabs(xOut) < 0.0000001) ? 0 : xOut;
0996     yOut = (fabs(yOut) < 0.0000001) ? 0 : yOut;
0997 
0998     // cout<<"Post rotation: "<<xOut<<" "<<yOut<<endl;
0999 
1000     return {xOut,yOut};
1001 }
1002 
1003 
1004 double RestDist::get_clu_eta(std::vector<double> vertex, std::vector<double> clu_pos)
1005 {
1006     double correct_x = clu_pos[0] - vertex[0];
1007     double correct_y = clu_pos[1] - vertex[1];
1008     double correct_z = clu_pos[2] - vertex[2];
1009     double clu_r = sqrt(pow(correct_x,2) + pow(correct_y,2));
1010     // cout<<"correct info : "<<correct_x<<" "<<correct_y<<" "<<correct_z<<" "<<clu_r<<endl;
1011 
1012     return -0.5 * TMath::Log((sqrt(pow(correct_z,2)+pow(clu_r,2))-(correct_z)) / (sqrt(pow(correct_z,2)+pow(clu_r,2))+(correct_z)));
1013 }
1014 
1015 double RestDist::calculateAngleBetweenVectors(double x1, double y1, double x2, double y2, double targetX, double targetY) {
1016     // Calculate the vectors vector_1 (point_1 to point_2) and vector_2 (point_1 to target)
1017     double vector1X = x2 - x1;
1018     double vector1Y = y2 - y1;
1019 
1020     double vector2X = targetX - x1;
1021     double vector2Y = targetY - y1;
1022 
1023     // Calculate the cross product of vector_1 and vector_2 (z-component)
1024     double crossProduct = vector1X * vector2Y - vector1Y * vector2X;
1025     
1026     // cout<<" crossProduct : "<<crossProduct<<endl;
1027 
1028     // Calculate the magnitudes of vector_1 and vector_2
1029     double magnitude1 = std::sqrt(vector1X * vector1X + vector1Y * vector1Y);
1030     double magnitude2 = std::sqrt(vector2X * vector2X + vector2Y * vector2Y);
1031 
1032     // Calculate the angle in radians using the inverse tangent of the cross product and dot product
1033     double dotProduct = vector1X * vector2X + vector1Y * vector2Y;
1034 
1035     double angleInRadians = std::atan2(std::abs(crossProduct), dotProduct);
1036     // Convert the angle from radians to degrees and return it
1037     double angleInDegrees = angleInRadians * 180.0 / M_PI;
1038     
1039     double angleInRadians_new = std::asin( crossProduct/(magnitude1*magnitude2) );
1040     double angleInDegrees_new = angleInRadians_new * 180.0 / M_PI;
1041     
1042     // cout<<"angle : "<<angleInDegrees_new<<endl;
1043 
1044     double DCA_distance = sin(angleInRadians_new) * magnitude2;
1045 
1046     return DCA_distance;
1047 }
[ 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!