Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "TrackletHistogramNew.h"
0002 
0003 TrackletHistogramNew::TrackletHistogramNew(
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     std::pair<bool, TH1D*> vtxZReweight_in,
0015     bool BcoFullDiffCut_in,
0016     bool INTT_vtxZ_QA_in,
0017     std::pair<bool, std::pair<double, double>> isClusQA_in,
0018     bool HaveGeoOffsetTag_in,
0019     std::pair<bool, int> SetRandomHits_in,
0020     bool RandInttZ_in,
0021     bool ColMulMask_in
0022 ) : ClusHistogram(
0023     process_id_in,
0024     runnumber_in,
0025     run_nEvents_in,
0026     input_directory_in,
0027     input_file_name_in,
0028     output_directory_in,
0029 
0030     output_file_name_suffix_in,
0031     vertexXYIncm_in,
0032 
0033     vtxZReweight_in,
0034     BcoFullDiffCut_in,
0035     INTT_vtxZ_QA_in,
0036     isClusQA_in,
0037     HaveGeoOffsetTag_in,
0038     SetRandomHits_in,
0039     RandInttZ_in,
0040     ColMulMask_in,
0041 
0042     1 // note : constructor type
0043 ){
0044     track_gr = new TGraphErrors();
0045     fit_rz = new TF1("fit_rz","pol1",-1000,1000);
0046 
0047     PrepareOutPutFileName();
0048     PrepareOutPutRootFile();
0049     PrepareHistograms();
0050 }
0051 
0052 void TrackletHistogramNew::PrepareOutPutFileName()
0053 {
0054     std::string job_index = std::to_string( process_id );
0055     int job_index_len = 5;
0056     job_index.insert(0, job_index_len - job_index.size(), '0');
0057 
0058     std::string runnumber_str = std::to_string( runnumber );
0059     if (runnumber != -1){
0060         int runnumber_str_len = 8;
0061         runnumber_str.insert(0, runnumber_str_len - runnumber_str.size(), '0');
0062     }
0063 
0064     if (output_file_name_suffix.size() > 0 && output_file_name_suffix[0] != '_') {
0065         output_file_name_suffix = "_" + output_file_name_suffix;
0066     }
0067 
0068     output_filename = "TrackHist";
0069     output_filename = (runnumber != -1) ? "Data_" + output_filename : "MC_" + output_filename;
0070     output_filename += (vtxZReweight.first) ? "_vtxZReweight" : "";
0071     output_filename += (BcoFullDiffCut && runnumber != -1) ? "_BcoFullDiffCut" : "";
0072     output_filename += (INTT_vtxZ_QA) ? "_VtxZQA" : "";
0073     output_filename += (isClusQA.first) ? Form("_ClusQAAdc%.0fPhiSize%.0f",isClusQA.second.first,isClusQA.second.second) : "";
0074     output_filename += (ColMulMask) ? "_ColMulMask" : "";
0075     output_filename += (HaveGeoOffsetTag) ? "_GeoOffset" : "";
0076     output_filename += (SetRandomHits.first) ? Form("_Rand%dHits",SetRandomHits.second) : "";
0077     output_filename += (RandInttZ) ? "_RandInttZ" : "";
0078     output_filename += output_file_name_suffix;
0079     output_filename += (runnumber != -1) ? Form("_%s_%s.root",runnumber_str.c_str(),job_index.c_str()) : Form("_%s.root",job_index.c_str());
0080 }
0081 
0082 void TrackletHistogramNew::PrepareOutPutRootFile()
0083 {
0084     file_out = new TFile(Form("%s/%s", output_directory.c_str(), output_filename.c_str()), "RECREATE");
0085     tree_out_par = new TTree("tree_par", "used parameters");
0086 
0087     tree_out_par -> Branch("process_id", &process_id);
0088     tree_out_par -> Branch("runnumber", &runnumber);
0089     tree_out_par -> Branch("vtxZReweight", &vtxZReweight.first);
0090     tree_out_par -> Branch("BcoFullDiffCut", &BcoFullDiffCut);
0091     tree_out_par -> Branch("INTT_vtxZ_QA", &INTT_vtxZ_QA);
0092     tree_out_par -> Branch("isClusQA", &isClusQA.first);
0093     tree_out_par -> Branch("isClusQA_adc", &isClusQA.second.first);
0094     tree_out_par -> Branch("isClusQA_phiSize", &isClusQA.second.second);
0095 
0096     tree_out_par -> Branch("centrality_edges", &centrality_edges);
0097     tree_out_par -> Branch("nCentrality_bin", &nCentrality_bin);
0098 
0099     tree_out_par -> Branch("CentralityFineEdge_min", &CentralityFineEdge_min);
0100     tree_out_par -> Branch("CentralityFineEdge_max", &CentralityFineEdge_max);
0101     tree_out_par -> Branch("nCentralityFineBin", &nCentralityFineBin);
0102 
0103     tree_out_par -> Branch("EtaEdge_min", &EtaEdge_min);
0104     tree_out_par -> Branch("EtaEdge_max", &EtaEdge_max);
0105     tree_out_par -> Branch("nEtaBin", &nEtaBin);
0106 
0107     tree_out_par -> Branch("VtxZEdge_min", &VtxZEdge_min);
0108     tree_out_par -> Branch("VtxZEdge_max", &VtxZEdge_max);
0109     tree_out_par -> Branch("nVtxZBin", &nVtxZBin);
0110 
0111     tree_out_par -> Branch("typeA_sensorZID", &typeA_sensorZID);
0112     tree_out_par -> Branch("cut_GlobalMBDvtxZ", &cut_GlobalMBDvtxZ);   
0113     tree_out_par -> Branch("cut_vtxZDiff", &cut_vtxZDiff);
0114     tree_out_par -> Branch("cut_TrapezoidalFitWidth", &cut_TrapezoidalFitWidth);
0115     tree_out_par -> Branch("cut_TrapezoidalFWHM", &cut_TrapezoidalFWHM);
0116     tree_out_par -> Branch("cut_INTTvtxZError", &cut_INTTvtxZError);
0117 
0118     tree_out_par -> Branch("cut_InttBcoFullDIff_next", &cut_InttBcoFullDIff_next);
0119 
0120     tree_out_par -> Branch("DeltaPhiEdge_min", &DeltaPhiEdge_min);
0121     tree_out_par -> Branch("DeltaPhiEdge_max", &DeltaPhiEdge_max);
0122     tree_out_par -> Branch("nDeltaPhiBin", &nDeltaPhiBin);
0123 
0124     tree_out_par -> Branch("cut_bestPair_DeltaPhi", &cut_bestPair_DeltaPhi);
0125     tree_out_par -> Branch("cut_GoodProtoTracklet_DeltaPhi", &cut_GoodProtoTracklet_DeltaPhi);
0126 }
0127 
0128 void TrackletHistogramNew::PrepareHistograms()
0129 {   
0130 
0131     h1D_GoodColMap_ZId = new TH1D("h1D_GoodColMap_ZId","h1D_GoodColMap_ZId;ClusZ [cm];Entries",nZbin, Zmin, Zmax);
0132 
0133     h1D_INTTvtxZ = new TH1D("h1D_INTTvtxZ", "h1D_INTTvtxZ;INTT vtxZ [cm];Entries", 60, -60, 60);
0134     h1D_centrality_bin = new TH1D("h1D_centrality_bin","h1D_centrality_bin;Centrality [%];Entries",nCentrality_bin,&centrality_edges[0]); // note : the 0-5%
0135     h1D_centrality_bin_weighted = new TH1D("h1D_centrality_bin_weighted","h1D_centrality_bin_weighted;Centrality [%];Entries",nCentrality_bin,&centrality_edges[0]); // note : the 0-5%
0136 
0137     h1D_INTTvtxZ_FineBin = new TH1D("h1D_INTTvtxZ_FineBin", "h1D_INTTvtxZ_FineBin;INTT vtxZ [cm];Entries", 120, -60, 60);
0138     h1D_INTTvtxZ_FineBin_NoVtxZWeight = new TH1D("h1D_INTTvtxZ_FineBin_NoVtxZWeight", "h1D_INTTvtxZ_FineBin_NoVtxZWeight;INTT vtxZ [cm];Entries", 120, -60, 60);
0139     h2D_INTTvtxZFineBin_CentralityBin = new TH2D("h2D_INTTvtxZFineBin_CentralityBin", "h2D_INTTvtxZFineBin_CentralityBin;INTT vtxZ [cm];Centrality [%]", 120, -60, 60, nCentrality_bin, &centrality_edges[0]);
0140 
0141     h1D_eta_template = new TH1D("h1D_eta_template", "h1D_eta_template", nEtaBin, EtaEdge_min, EtaEdge_max); // note : coarse
0142     h1D_vtxz_template = new TH1D("h1D_vtxz_template", "h1D_vtxz_template", nVtxZBin, VtxZEdge_min, VtxZEdge_max); // note : coarse
0143 
0144     // note : for inclusive 
0145     h1D_PairDeltaEta_inclusive = new TH1D("h1D_PairDeltaEta_inclusive", "h1D_PairDeltaEta_inclusive;Pair #Delta#eta;Entries", nDeltaEtaBin, DeltaEtaEdge_min, DeltaEtaEdge_max);
0146     h1D_PairDeltaPhi_inclusive = new TH1D("h1D_PairDeltaPhi_inclusive", "h1D_PairDeltaPhi_inclusive;Pair #Delta#phi;Entries", nDeltaPhiBin, DeltaPhiEdge_min, DeltaPhiEdge_max);
0147 
0148     // note : cluster Eta
0149     h1D_Inner_ClusEta_INTTz = new TH1D("h1D_Inner_ClusEta_INTTz","h1D_Inner_ClusEta_INTTz;ClusEta (Inner, w.r.t INTTz);Entries", nEtaBin, EtaEdge_min, EtaEdge_max);
0150     h1D_Outer_ClusEta_INTTz = new TH1D("h1D_Outer_ClusEta_INTTz","h1D_Outer_ClusEta_INTTz;ClusEta (Outer, w.r.t INTTz);Entries", nEtaBin, EtaEdge_min, EtaEdge_max);
0151 
0152     std::vector<int> insert_check; insert_check.clear();
0153     bool isInserted = false;
0154 
0155     // note : for 1D
0156     h1D_map.clear();
0157 
0158     isInserted = h1D_map.insert( std::make_pair(
0159             "h1D_centrality",
0160             new TH1D("h1D_centrality", "h1D_centrality;Centrality [%];Entries", nCentralityFineBin, CentralityFineEdge_min, CentralityFineEdge_max)
0161         )
0162     ).second; insert_check.push_back(isInserted);
0163 
0164     for (int vtxz_bin = 0; vtxz_bin < nVtxZBin; vtxz_bin++)
0165     {
0166         isInserted = h1D_map.insert( std::make_pair(
0167                 Form("h1D_centrality_InttVtxZ%d", vtxz_bin),
0168                 new TH1D(Form("h1D_centrality_InttVtxZ%d", vtxz_bin), Form("h1D_centrality_InttVtxZ%d;Centrality [%];Entries",vtxz_bin), nCentralityFineBin, CentralityFineEdge_min, CentralityFineEdge_max)
0169             )
0170         ).second; insert_check.push_back(isInserted);
0171 
0172         // isInserted = h1D_map.insert( std::make_pair(
0173         //         Form("h1D_centrality_MBDVtxZ%d", vtxz_bin),
0174         //         new TH1D(Form("h1D_centrality_MBDVtxZ%d", vtxz_bin), Form("h1D_centrality_MBDVtxZ%d;Centrality [%];Entries",vtxz_bin), nCentralityFineBin, CentralityFineEdge_min, CentralityFineEdge_max)
0175         //     )
0176         // ).second; insert_check.push_back(isInserted);
0177     }
0178 
0179     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0180 
0181     for (int eta_bin = 0; eta_bin < nEtaBin; eta_bin++)
0182     {
0183         for (int vtxz_bin = 0; vtxz_bin < nVtxZBin; vtxz_bin++)
0184         {
0185             for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++)
0186             {
0187                 isInserted = h1D_map.insert( std::make_pair(
0188                         Form("h1D_DeltaPhi_Mbin%d_Eta%d_VtxZ%d", Mbin, eta_bin, vtxz_bin),
0189                         new TH1D(Form("h1D_DeltaPhi_Mbin%d_Eta%d_VtxZ%d", Mbin, eta_bin, vtxz_bin), Form("h1D_DeltaPhi_Mbin%d_Eta%d_VtxZ%d;Pair #Delta#phi [radian];Entries", Mbin, eta_bin, vtxz_bin), nDeltaPhiBin, DeltaPhiEdge_min, DeltaPhiEdge_max)
0190                     )
0191                 ).second; insert_check.push_back(isInserted);
0192 
0193                 isInserted = h1D_map.insert( std::make_pair(
0194                         Form("h1D_DeltaPhi_Mbin%d_Eta%d_VtxZ%d_rotated", Mbin, eta_bin, vtxz_bin),
0195                         new TH1D(Form("h1D_DeltaPhi_Mbin%d_Eta%d_VtxZ%d_rotated", Mbin, eta_bin, vtxz_bin), Form("h1D_DeltaPhi_Mbin%d_Eta%d_VtxZ%d_rotated;Pair #Delta#phi [radian];Entries", Mbin, eta_bin, vtxz_bin), nDeltaPhiBin, DeltaPhiEdge_min, DeltaPhiEdge_max)
0196                     )
0197                 ).second; insert_check.push_back(isInserted);
0198                 
0199 
0200 
0201                 isInserted = h1D_map.insert( std::make_pair(
0202                         Form("h1D_typeA_DeltaPhi_Mbin%d_Eta%d_VtxZ%d", Mbin, eta_bin, vtxz_bin),
0203                         new TH1D(Form("h1D_typeA_DeltaPhi_Mbin%d_Eta%d_VtxZ%d", Mbin, eta_bin, vtxz_bin), Form("h1D_typeA_DeltaPhi_Mbin%d_Eta%d_VtxZ%d;Pair #Delta#phi [radian] (type A);Entries", Mbin, eta_bin, vtxz_bin), nDeltaPhiBin, DeltaPhiEdge_min, DeltaPhiEdge_max)
0204                     )
0205                 ).second; insert_check.push_back(isInserted);
0206 
0207                 isInserted = h1D_map.insert( std::make_pair(
0208                         Form("h1D_typeA_DeltaPhi_Mbin%d_Eta%d_VtxZ%d_rotated", Mbin, eta_bin, vtxz_bin),
0209                         new TH1D(Form("h1D_typeA_DeltaPhi_Mbin%d_Eta%d_VtxZ%d_rotated", Mbin, eta_bin, vtxz_bin), Form("h1D_typeA_DeltaPhi_Mbin%d_Eta%d_VtxZ%d_rotated;Pair #Delta#phi [radian] (type A);Entries", Mbin, eta_bin, vtxz_bin), nDeltaPhiBin, DeltaPhiEdge_min, DeltaPhiEdge_max)
0210                     )
0211                 ).second; insert_check.push_back(isInserted);
0212                 
0213 
0214 
0215                 // note : for truth eta
0216                 if (runnumber == -1 && eta_bin == 0){
0217                     isInserted = h1D_map.insert( std::make_pair(
0218                             Form("h1D_TrueEta_Mbin%d_VtxZ%d", Mbin, vtxz_bin),
0219                             new TH1D(Form("h1D_TrueEta_Mbin%d_VtxZ%d", Mbin, vtxz_bin), Form("h1D_TrueEta_Mbin%d_VtxZ%d;PHG4Particle #eta;Entries", Mbin, vtxz_bin), nEtaBin, EtaEdge_min, EtaEdge_max)
0220                         )
0221                     ).second; insert_check.push_back(isInserted);
0222                 }
0223             } // note : end Mbin
0224         }
0225     }// note : end of Delta_Phi hist
0226 
0227     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0228 
0229     // note : for the best pairs
0230     for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++)
0231     {
0232         for (int vtxz_bin = 0; vtxz_bin < nVtxZBin; vtxz_bin++){
0233 
0234             isInserted = h1D_map.insert( std::make_pair(
0235                     Form("h1D_BestPair_DeltaEta_Mbin%d_VtxZ%d",Mbin, vtxz_bin),
0236                     new TH1D(Form("h1D_BestPair_DeltaEta_Mbin%d_VtxZ%d",Mbin, vtxz_bin), Form("h1D_BestPair_DeltaEta_Mbin%d_VtxZ%d;Pair #Delta#eta;Entries",Mbin, vtxz_bin), nDeltaEtaBin, DeltaEtaEdge_min, DeltaEtaEdge_max)
0237                 )
0238             ).second; insert_check.push_back(isInserted);
0239 
0240             isInserted = h1D_map.insert( std::make_pair(
0241                     Form("h1D_BestPair_DeltaPhi_Mbin%d_VtxZ%d",Mbin, vtxz_bin),
0242                     new TH1D(Form("h1D_BestPair_DeltaPhi_Mbin%d_VtxZ%d",Mbin, vtxz_bin), Form("h1D_BestPair_DeltaPhi_Mbin%d_VtxZ%d;Pair #Delta#phi [radian];Entries",Mbin, vtxz_bin), nDeltaPhiBin, DeltaPhiEdge_min, DeltaPhiEdge_max)
0243                 )
0244             ).second; insert_check.push_back(isInserted);
0245 
0246             isInserted = h1D_map.insert( std::make_pair(
0247                     Form("h1D_BestPair_ClusPhiSize_Mbin%d_VtxZ%d",Mbin, vtxz_bin),
0248                     new TH1D(Form("h1D_BestPair_ClusPhiSize_Mbin%d_VtxZ%d",Mbin, vtxz_bin), Form("h1D_BestPair_ClusPhiSize_Mbin%d_VtxZ%d;ClusPhiSize;Entries",Mbin, vtxz_bin), 128,0,128)
0249                 )
0250             ).second; insert_check.push_back(isInserted);
0251 
0252             isInserted = h1D_map.insert( std::make_pair(
0253                     Form("h1D_BestPair_ClusAdc_Mbin%d_VtxZ%d",Mbin, vtxz_bin),
0254                     new TH1D(Form("h1D_BestPair_ClusAdc_Mbin%d_VtxZ%d",Mbin, vtxz_bin), Form("h1D_BestPair_ClusAdc_Mbin%d_VtxZ%d;ClusAdc;Entries",Mbin, vtxz_bin), 200,0,18000)
0255                 )
0256             ).second; insert_check.push_back(isInserted);
0257 
0258             
0259             // note : typeA : ------------------------------------------------------------------------------------------------------------------------------------------------------------------
0260             isInserted = h1D_map.insert( std::make_pair(
0261                     Form("h1D_typeA_BestPair_DeltaEta_Mbin%d_VtxZ%d",Mbin, vtxz_bin),
0262                     new TH1D(Form("h1D_typeA_BestPair_DeltaEta_Mbin%d_VtxZ%d",Mbin, vtxz_bin), Form("h1D_typeA_BestPair_DeltaEta_Mbin%d_VtxZ%d;Pair #Delta#eta (type A);Entries",Mbin, vtxz_bin), nDeltaEtaBin, DeltaEtaEdge_min, DeltaEtaEdge_max)
0263                 )
0264             ).second; insert_check.push_back(isInserted);
0265 
0266             isInserted = h1D_map.insert( std::make_pair(
0267                     Form("h1D_typeA_BestPair_DeltaPhi_Mbin%d_VtxZ%d",Mbin, vtxz_bin),
0268                     new TH1D(Form("h1D_typeA_BestPair_DeltaPhi_Mbin%d_VtxZ%d",Mbin, vtxz_bin), Form("h1D_typeA_BestPair_DeltaPhi_Mbin%d_VtxZ%d;Pair #Delta#phi [radian] (type A);Entries",Mbin, vtxz_bin), nDeltaPhiBin, DeltaPhiEdge_min, DeltaPhiEdge_max)
0269                 )
0270             ).second; insert_check.push_back(isInserted);
0271 
0272             isInserted = h1D_map.insert( std::make_pair(
0273                     Form("h1D_typeA_BestPair_ClusPhiSize_Mbin%d_VtxZ%d",Mbin, vtxz_bin),
0274                     new TH1D(Form("h1D_typeA_BestPair_ClusPhiSize_Mbin%d_VtxZ%d",Mbin, vtxz_bin), Form("h1D_typeA_BestPair_ClusPhiSize_Mbin%d_VtxZ%d;ClusPhiSize (type A);Entries",Mbin, vtxz_bin), 128,0,128)
0275                 )
0276             ).second; insert_check.push_back(isInserted);
0277 
0278             isInserted = h1D_map.insert( std::make_pair(
0279                     Form("h1D_typeA_BestPair_ClusAdc_Mbin%d_VtxZ%d",Mbin, vtxz_bin),
0280                     new TH1D(Form("h1D_typeA_BestPair_ClusAdc_Mbin%d_VtxZ%d",Mbin, vtxz_bin), Form("h1D_typeA_BestPair_ClusAdc_Mbin%d_VtxZ%d;ClusAdc (type A);Entries",Mbin, vtxz_bin), 200,0,18000)
0281                 )
0282             ).second; insert_check.push_back(isInserted);
0283         }
0284 
0285     }
0286 
0287 
0288     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0289 
0290     // note : for 2D
0291     h2D_map.clear();
0292 
0293     // note : truth nPHG4Particle in Eta-VtxZ
0294     if (runnumber == -1){
0295         for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++)
0296         {
0297             isInserted = h2D_map.insert( std::make_pair(
0298                     Form("h2D_TrueEtaVtxZ_Mbin%d", Mbin),
0299                     new TH2D(Form("h2D_TrueEtaVtxZ_Mbin%d", Mbin), Form("h2D_TrueEtaVtxZ_Mbin%d;PHG4Particle #eta;INTT vtxZ [cm]", Mbin), nEtaBin, EtaEdge_min, EtaEdge_max, nVtxZBin, VtxZEdge_min, VtxZEdge_max)
0300                 )
0301             ).second; insert_check.push_back(isInserted);
0302 
0303             isInserted = h2D_map.insert( std::make_pair(
0304                     Form("h2D_TrueEtaVtxZ_Mbin%d_FineBin", Mbin),
0305                     new TH2D(Form("h2D_TrueEtaVtxZ_Mbin%d_FineBin", Mbin), Form("h2D_TrueEtaVtxZ_Mbin%d_FineBin;PHG4Particle #eta;INTT vtxZ [cm]", Mbin), 540, EtaEdge_min, EtaEdge_max, nVtxZBin_FineBin, VtxZEdge_min, VtxZEdge_max)
0306                 )
0307             ).second; insert_check.push_back(isInserted);
0308         }
0309 
0310         isInserted = h1D_map.insert( std::make_pair(
0311                 "debug_h1D_NHadron_Inclusive100",
0312                 new TH1D("debug_h1D_NHadron_Inclusive100", "debug_h1D_NHadron_Inclusive100;N Hadrons;Entries", 200, 0, 10000)
0313             )
0314         ).second; insert_check.push_back(isInserted);
0315 
0316         isInserted = h1D_map.insert( std::make_pair(
0317                 "debug_h1D_NHadron_OneEtaBin_Inclusive100",
0318                 new TH1D("debug_h1D_NHadron_OneEtaBin_Inclusive100", "debug_h1D_NHadron_OneEtaBin_Inclusive100;N Hadrons;Entries", 200, 0, 400)
0319             )
0320         ).second; insert_check.push_back(isInserted);
0321 
0322         isInserted = h2D_map.insert( std::make_pair(
0323                 Form("h2D_TrueEvtCount_vtxZCentrality"),
0324                 new TH2D(Form("h2D_TrueEvtCount_vtxZCentrality"), Form("h2D_TrueEvtCount_vtxZCentrality;INTT vtxZ [cm];Centrality [%]"), nVtxZBin, VtxZEdge_min, VtxZEdge_max, nCentrality_bin, &centrality_edges[0]) // note : 0 - 5%
0325             )
0326         ).second; insert_check.push_back(isInserted);
0327     }
0328 
0329     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0330 
0331     // note : for good proto tracklet in Eta-VtxZ
0332     for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
0333         
0334         isInserted = h2D_map.insert( std::make_pair(
0335                 Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d", Mbin),
0336                 new TH2D(Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d", Mbin), Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d;Pair #eta;INTT vtxZ [cm]", Mbin), nEtaBin, EtaEdge_min, EtaEdge_max, nVtxZBin, VtxZEdge_min, VtxZEdge_max) 
0337             )
0338         ).second; insert_check.push_back(isInserted);
0339 
0340         isInserted = h2D_map.insert( std::make_pair(
0341                 Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d", Mbin),
0342                 new TH2D(Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d", Mbin), Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d;Pair #eta (type A);INTT vtxZ [cm]", Mbin), nEtaBin, EtaEdge_min, EtaEdge_max, nVtxZBin, VtxZEdge_min, VtxZEdge_max) 
0343             )
0344         ).second; insert_check.push_back(isInserted);
0345 
0346         // note : normal fine bin
0347         isInserted = h2D_map.insert( std::make_pair(
0348                 Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin", Mbin),
0349                 new TH2D(Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin", Mbin), Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin;Pair #eta;INTT vtxZ [cm]", Mbin), 540, EtaEdge_min, EtaEdge_max, nVtxZBin_FineBin, VtxZEdge_min, VtxZEdge_max) 
0350             )
0351         ).second; insert_check.push_back(isInserted);
0352 
0353         isInserted = h2D_map.insert( std::make_pair(
0354                 Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin", Mbin),
0355                 new TH2D(Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin", Mbin), Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin;Pair #eta (type A);INTT vtxZ [cm]", Mbin), 540, EtaEdge_min, EtaEdge_max, nVtxZBin_FineBin, VtxZEdge_min, VtxZEdge_max) 
0356             )
0357         ).second; insert_check.push_back(isInserted);
0358         
0359 
0360         // note : rotated
0361         isInserted = h2D_map.insert( std::make_pair(
0362                 Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated", Mbin),
0363                 new TH2D(Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated", Mbin), Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated;Pair #eta;INTT vtxZ [cm]", Mbin), nEtaBin, EtaEdge_min, EtaEdge_max, nVtxZBin, VtxZEdge_min, VtxZEdge_max) 
0364             )
0365         ).second; insert_check.push_back(isInserted);
0366 
0367         isInserted = h2D_map.insert( std::make_pair(
0368                 Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated", Mbin),
0369                 new TH2D(Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated", Mbin), Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated;Pair #eta (type A);INTT vtxZ [cm]", Mbin), nEtaBin, EtaEdge_min, EtaEdge_max, nVtxZBin, VtxZEdge_min, VtxZEdge_max) 
0370             )
0371         ).second; insert_check.push_back(isInserted);
0372 
0373         // note : FineBin
0374         isInserted = h2D_map.insert( std::make_pair(
0375                 Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated_FineBin", Mbin),
0376                 new TH2D(Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated_FineBin", Mbin), Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated_FineBin;Pair #eta;INTT vtxZ [cm]", Mbin), 540, EtaEdge_min, EtaEdge_max, nVtxZBin_FineBin, VtxZEdge_min, VtxZEdge_max) 
0377             )
0378         ).second; insert_check.push_back(isInserted);
0379 
0380         isInserted = h2D_map.insert( std::make_pair(
0381                 Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated_FineBin", Mbin),
0382                 new TH2D(Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated_FineBin", Mbin), Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated_FineBin;Pair #eta (type A);INTT vtxZ [cm]", Mbin), 540, EtaEdge_min, EtaEdge_max, nVtxZBin_FineBin, VtxZEdge_min, VtxZEdge_max) 
0383             )
0384         ).second; insert_check.push_back(isInserted);
0385     }
0386 
0387 
0388     isInserted = h2D_map.insert( std::make_pair(
0389             Form("h2D_RecoEvtCount_vtxZCentrality"),
0390             new TH2D(Form("h2D_RecoEvtCount_vtxZCentrality"), Form("h2D_RecoEvtCount_vtxZCentrality;INTT vtxZ [cm];Centrality [%]"), nVtxZBin, VtxZEdge_min, VtxZEdge_max, nCentrality_bin, &centrality_edges[0]) // note : 0 - 5%
0391         )
0392     ).second; insert_check.push_back(isInserted);
0393     
0394     // note : best pair in Eta-VtxZ
0395     for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
0396         isInserted = h2D_map.insert( std::make_pair(
0397                 Form("h2D_BestPairEtaVtxZ_Mbin%d",Mbin),
0398                 new TH2D(Form("h2D_BestPairEtaVtxZ_Mbin%d",Mbin), Form("h2D_BestPairEtaVtxZ_Mbin%d;Tracklet #eta;INTT vtxZ [cm]",Mbin), nEtaBin, EtaEdge_min, EtaEdge_max, nVtxZBin, VtxZEdge_min, VtxZEdge_max)
0399             )
0400         ).second; insert_check.push_back(isInserted);
0401 
0402         isInserted = h2D_map.insert( std::make_pair(
0403                 Form("h2D_BestPairEtaVtxZ_Mbin%d_FineBin",Mbin),
0404                 new TH2D(Form("h2D_BestPairEtaVtxZ_Mbin%d_FineBin",Mbin), Form("h2D_BestPairEtaVtxZ_Mbin%d_FineBin;Tracklet #eta;INTT vtxZ [cm]",Mbin), 540, EtaEdge_min, EtaEdge_max, nVtxZBin_FineBin, VtxZEdge_min, VtxZEdge_max)
0405             )
0406         ).second; insert_check.push_back(isInserted);
0407 
0408         isInserted = h2D_map.insert( std::make_pair(
0409                 Form("h2D_typeA_BestPairEtaVtxZ_Mbin%d",Mbin),
0410                 new TH2D(Form("h2D_typeA_BestPairEtaVtxZ_Mbin%d",Mbin), Form("h2D_typeA_BestPairEtaVtxZ_Mbin%d;Tracklet #eta (type A);INTT vtxZ [cm]",Mbin), nEtaBin, EtaEdge_min, EtaEdge_max, nVtxZBin, VtxZEdge_min, VtxZEdge_max)
0411             )
0412         ).second; insert_check.push_back(isInserted);
0413 
0414         isInserted = h2D_map.insert( std::make_pair(
0415                 Form("h2D_typeA_BestPairEtaVtxZ_Mbin%d_FineBin",Mbin),
0416                 new TH2D(Form("h2D_typeA_BestPairEtaVtxZ_Mbin%d_FineBin",Mbin), Form("h2D_typeA_BestPairEtaVtxZ_Mbin%d_FineBin;Tracklet #eta (type A);INTT vtxZ [cm]",Mbin), 540, EtaEdge_min, EtaEdge_max, nVtxZBin_FineBin, VtxZEdge_min, VtxZEdge_max)
0417             )
0418         ).second; insert_check.push_back(isInserted);
0419     }
0420 
0421     
0422 
0423 
0424     // note : vtxZ - centrality
0425     isInserted = h2D_map.insert( std::make_pair(
0426             Form("h2D_InttVtxZ_Centrality"),
0427             new TH2D(Form("h2D_InttVtxZ_Centrality"), Form("h2D_InttVtxZ_Centrality;INTT vtxZ [cm];Centrality"), nVtxZBin, VtxZEdge_min, VtxZEdge_max, nCentralityFineBin, CentralityFineEdge_min, CentralityFineEdge_max)
0428         )
0429     ).second; insert_check.push_back(isInserted);
0430 
0431     // isInserted = h2D_map.insert( std::make_pair(
0432     //         Form("h2D_MBDVtxZ_Centrality"),
0433     //         new TH2D(Form("h2D_MBDVtxZ_Centrality"), Form("h2D_MBDVtxZ_Centrality;INTT vtxZ [cm];Centrality"), nVtxZBin, VtxZEdge_min, VtxZEdge_max, nCentralityFineBin, CentralityFineEdge_min, CentralityFineEdge_max)
0434     //     )
0435     // ).second; insert_check.push_back(isInserted);
0436     
0437 
0438     for (int isInserted_tag : insert_check){
0439         if (isInserted_tag == 0){
0440             std::cout<<"Histogram insertion failed"<<std::endl;
0441             exit(1);
0442         }
0443     }
0444 }
0445 
0446 void TrackletHistogramNew::EvtCleanUp()
0447 {
0448     Used_Clus_index_vec.clear();
0449     Pair_DeltaPhi_vec.clear();
0450 
0451     evt_sPH_inner_nocolumn_vec.clear();
0452     evt_sPH_outer_nocolumn_vec.clear();
0453 
0454     evt_TrackletPair_vec.clear();
0455     evt_TrackletPairRotate_vec.clear();
0456 
0457     inner_clu_phi_map.clear();
0458     outer_clu_phi_map.clear();
0459     inner_clu_phi_map = std::vector<std::vector<std::pair<bool,ClusHistogram::clu_info>>>(360);
0460     outer_clu_phi_map = std::vector<std::vector<std::pair<bool,ClusHistogram::clu_info>>>(360);
0461 }
0462 
0463 void TrackletHistogramNew::GetTrackletPair(std::vector<pair_str> &input_TrackletPair_vec, bool isRotated)
0464 {
0465     input_TrackletPair_vec.clear();
0466 
0467     inner_clu_phi_map.clear();
0468     outer_clu_phi_map.clear();
0469     inner_clu_phi_map = std::vector<std::vector<std::pair<bool,ClusHistogram::clu_info>>>(360);
0470     outer_clu_phi_map = std::vector<std::vector<std::pair<bool,ClusHistogram::clu_info>>>(360);
0471 
0472     if (INTTvtxZ != INTTvtxZ || INTTvtxZError != INTTvtxZError) {return;}
0473 
0474     std::vector<ClusHistogram::clu_info> temp_evt_sPH_inner_nocolumn_vec = (isRotated) ? GetRotatedClusterVec(evt_sPH_inner_nocolumn_vec) : evt_sPH_inner_nocolumn_vec;
0475 
0476     for (int inner_i = 0; inner_i < int(temp_evt_sPH_inner_nocolumn_vec.size()); inner_i++) {
0477       double Clus_InnerPhi_Offset = (temp_evt_sPH_inner_nocolumn_vec[inner_i].y - vertexXYIncm.second < 0) ? atan2(temp_evt_sPH_inner_nocolumn_vec[inner_i].y - vertexXYIncm.second, temp_evt_sPH_inner_nocolumn_vec[inner_i].x - vertexXYIncm.first) * (180./TMath::Pi()) + 360 : atan2(temp_evt_sPH_inner_nocolumn_vec[inner_i].y - vertexXYIncm.second, temp_evt_sPH_inner_nocolumn_vec[inner_i].x - vertexXYIncm.first) * (180./TMath::Pi());
0478       inner_clu_phi_map[ int(Clus_InnerPhi_Offset) ].push_back({false,temp_evt_sPH_inner_nocolumn_vec[inner_i]});
0479     }
0480     for (int outer_i = 0; outer_i < int(evt_sPH_outer_nocolumn_vec.size()); outer_i++) {
0481         double Clus_OuterPhi_Offset = (evt_sPH_outer_nocolumn_vec[outer_i].y - vertexXYIncm.second < 0) ? atan2(evt_sPH_outer_nocolumn_vec[outer_i].y - vertexXYIncm.second, evt_sPH_outer_nocolumn_vec[outer_i].x - vertexXYIncm.first) * (180./TMath::Pi()) + 360 : atan2(evt_sPH_outer_nocolumn_vec[outer_i].y - vertexXYIncm.second, evt_sPH_outer_nocolumn_vec[outer_i].x - vertexXYIncm.first) * (180./TMath::Pi());
0482         outer_clu_phi_map[ int(Clus_OuterPhi_Offset) ].push_back({false,evt_sPH_outer_nocolumn_vec[outer_i]});
0483     }
0484 
0485 
0486     for (int inner_phi_i = 0; inner_phi_i < 360; inner_phi_i++) // note : each phi cell (1 degree)
0487     {
0488         // note : N cluster in this phi cell
0489         for (int inner_phi_clu_i = 0; inner_phi_clu_i < inner_clu_phi_map[inner_phi_i].size(); inner_phi_clu_i++)
0490         {
0491             if (inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].first == true) {continue;}
0492 
0493             ClusHistogram::clu_info inner_clu = inner_clu_phi_map[inner_phi_i][inner_phi_clu_i].second;
0494 
0495             double Clus_InnerPhi_Offset_radian = atan2(inner_clu.y - vertexXYIncm.second, inner_clu.x - vertexXYIncm.first);
0496             double Clus_InnerPhi_Offset = (inner_clu.y - vertexXYIncm.second < 0) ? Clus_InnerPhi_Offset_radian * (180./TMath::Pi()) + 360 : Clus_InnerPhi_Offset_radian * (180./TMath::Pi());
0497 
0498             // todo: change the outer phi scan range
0499             // note : the outer phi index, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5
0500             for (int scan_i = -5; scan_i < 6; scan_i++)
0501             {
0502                 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;
0503 
0504                 // note : N clusters in that outer phi cell
0505                 for (int outer_phi_clu_i = 0; outer_phi_clu_i < outer_clu_phi_map[true_scan_i].size(); outer_phi_clu_i++)
0506                 {
0507                     if (outer_clu_phi_map[true_scan_i][outer_phi_clu_i].first == true) {continue;}
0508 
0509                     ClusHistogram::clu_info outer_clu = outer_clu_phi_map[true_scan_i][outer_phi_clu_i].second;
0510 
0511                     double Clus_OuterPhi_Offset_radian = atan2(outer_clu.y - vertexXYIncm.second, outer_clu.x - vertexXYIncm.first);
0512                     double Clus_OuterPhi_Offset = (outer_clu.y - vertexXYIncm.second < 0) ? Clus_OuterPhi_Offset_radian * (180./TMath::Pi()) + 360 : Clus_OuterPhi_Offset_radian * (180./TMath::Pi());
0513                     double delta_phi = get_delta_phi(Clus_InnerPhi_Offset, Clus_OuterPhi_Offset);
0514                     // note : degree to radian 
0515                     double delta_phi_radian = delta_phi * (TMath::Pi() / 180.);
0516                     
0517                     // if (fabs(delta_phi) > 5.72) {continue;}
0518                     // if (fabs(delta_phi) > 3.5) {continue;}
0519 
0520                     double inner_clu_eta = get_clu_eta({vertexXYIncm.first, vertexXYIncm.second, INTTvtxZ},{inner_clu.x, inner_clu.y, inner_clu.z});
0521                     double outer_clu_eta = get_clu_eta({vertexXYIncm.first, vertexXYIncm.second, INTTvtxZ},{outer_clu.x, outer_clu.y, outer_clu.z});
0522                     double delta_eta = inner_clu_eta - outer_clu_eta;
0523 
0524                     std::pair<double,double> z_range_info = Get_possible_zvtx( 
0525                         0., // get_radius(vertexXYIncm.first,vertexXYIncm.second), 
0526                         {
0527                             get_radius(inner_clu.x - vertexXYIncm.first, inner_clu.y - vertexXYIncm.second), 
0528                             inner_clu.z,
0529                             double(inner_clu.sensorZID)
0530                         }, // note : unsign radius
0531                         
0532                         {
0533                             get_radius(outer_clu.x - vertexXYIncm.first, outer_clu.y - vertexXYIncm.second), 
0534                             outer_clu.z,
0535                             double(outer_clu.sensorZID)
0536                         }  // note : unsign radius
0537                     );
0538 
0539                     // note : this is a cut to constraint on the z vertex, only if the tracklets with the range that covers the z vertex can pass this cut
0540                     if (z_range_info.first - z_range_info.second > INTTvtxZ + INTTvtxZError || z_range_info.first + z_range_info.second < INTTvtxZ - INTTvtxZError) {continue;}
0541 
0542                     // note : do the fill here (find the best match outer cluster with the inner cluster )
0543                     // std::pair<double,double> Get_eta_pair = Get_eta(
0544                     //     {0., INTTvtxZ,INTTvtxZError},
0545                     //     {
0546                     //         get_radius(inner_clu.x - vertexXYIncm.first, inner_clu.y - vertexXYIncm.second), 
0547                     //         inner_clu.z,
0548                     //         double(inner_clu.sensorZID)
0549                     //     },
0550                     //     {
0551                     //         get_radius(outer_clu.x - vertexXYIncm.first, outer_clu.y - vertexXYIncm.second), 
0552                     //         outer_clu.z,
0553                     //         double(outer_clu.sensorZID)
0554                     //     }
0555                     // );  
0556 
0557                     // if (isRotated == false && fabs(delta_phi_radian) < 0.021 && INTTvtxZ > 30 && Get_eta_pair.second > 0){
0558                     //     std::cout<<"Check, vertex: {"<<vertexXYIncm.first<<", "<<vertexXYIncm.second<<", "<<INTTvtxZ<<"}, INTTvtxZError: "<<INTTvtxZError<<std::endl;
0559                     //     std::cout<<"Inner cluster Pos: {"<<inner_clu.x<<", "<<inner_clu.y<<", "<<inner_clu.z<<"}, eta_INTTz: "<<inner_clu_eta<<", sensorID: "<<inner_clu.sensorZID<<std::endl;
0560                     //     std::cout<<"Outer cluster Pos: {"<<outer_clu.x<<", "<<outer_clu.y<<", "<<outer_clu.z<<"}, eta_INTTz: "<<outer_clu_eta<<", sensorID: "<<outer_clu.sensorZID<<std::endl;
0561                     //     std::cout<<"DeltaPhi: "<<delta_phi_radian<<", DeltaEta: "<<delta_eta<<", PairEta: "<<(inner_clu_eta + outer_clu_eta) / 2.<<", PairEtaFit: "<<Get_eta_pair.second<<std::endl;
0562                     //     std::cout<<std::endl;
0563                     // }
0564 
0565                     pair_str temp_pair_str;
0566                     temp_pair_str.delta_phi = delta_phi_radian;
0567                     temp_pair_str.delta_eta = delta_eta;
0568                     temp_pair_str.pair_eta_num = (inner_clu_eta + outer_clu_eta) / 2.;
0569                     temp_pair_str.pair_eta_fit = std::nan("");//Get_eta_pair.second; // todo: Cancel the eta fitting
0570 
0571                     // temp_pair_str.inner_phi_id = inner_clu.ladderPhiID;
0572                     // temp_pair_str.inner_layer_id = inner_clu.layerID;
0573                     temp_pair_str.inner_zid = inner_clu.sensorZID;
0574                     temp_pair_str.inner_phi_size = inner_clu.phi_size;
0575                     temp_pair_str.inner_adc = inner_clu.adc;
0576                     temp_pair_str.inner_x = inner_clu.x;
0577                     temp_pair_str.inner_y = inner_clu.y;
0578                     // temp_pair_str.inner_z = inner_clu.z;
0579                     // temp_pair_str.inner_phi = Clus_InnerPhi_Offset_radian;
0580                     // temp_pair_str.inner_eta = inner_clu_eta;
0581                     temp_pair_str.inner_index = inner_clu.index;
0582                     
0583                     // temp_pair_str.outer_phi_id = outer_clu.ladderPhiID;
0584                     // temp_pair_str.outer_layer_id = outer_clu.layerID;
0585                     temp_pair_str.outer_zid = outer_clu.sensorZID;
0586                     temp_pair_str.outer_phi_size = outer_clu.phi_size;
0587                     temp_pair_str.outer_adc = outer_clu.adc;
0588                     temp_pair_str.outer_x = outer_clu.x;
0589                     temp_pair_str.outer_y = outer_clu.y;
0590                     // temp_pair_str.outer_z = outer_clu.z;
0591                     // temp_pair_str.outer_phi = Clus_OuterPhi_Offset_radian;
0592                     // temp_pair_str.outer_eta = outer_clu_eta;
0593                     temp_pair_str.outer_index = outer_clu.index;
0594 
0595                     input_TrackletPair_vec.push_back(temp_pair_str);
0596 
0597                 }
0598             }
0599         }
0600     }
0601 
0602 
0603 }
0604 
0605 void TrackletHistogramNew::FillPairs(std::vector<pair_str> input_TrackletPair_vec, bool isRotated, int Mbin_in, int vtxz_bin_in, double vtxZ_weight_in, int eID_in)
0606 {
0607     Pair_DeltaPhi_vec.clear();
0608     Used_Clus_index_vec.clear();
0609 
0610     std::string rotated_text = (isRotated) ? "_rotated" : "";
0611 
0612     for (int pair_i = 0; pair_i < input_TrackletPair_vec.size(); pair_i++) // note : either non-rotated or rotated
0613     {
0614         pair_str this_pair = input_TrackletPair_vec[pair_i];
0615 
0616         if (isRotated == false) {Pair_DeltaPhi_vec.push_back(fabs(this_pair.delta_phi));}
0617 
0618         int eta_bin = h1D_eta_template -> Fill(this_pair.pair_eta_num);
0619         eta_bin = (eta_bin == -1) ? -1 : eta_bin - 1;
0620 
0621         if (eta_bin == -1) {continue;}
0622 
0623         if (isRotated == false){
0624             h1D_PairDeltaEta_inclusive -> Fill(this_pair.delta_eta, vtxZ_weight_in);
0625             h1D_PairDeltaPhi_inclusive -> Fill(this_pair.delta_phi, vtxZ_weight_in);
0626         }
0627 
0628         // note : for both rotated and non-rotated
0629         h1D_map[Form("h1D_DeltaPhi_Mbin%d_Eta%d_VtxZ%d%s", Mbin_in, eta_bin, vtxz_bin_in, rotated_text.c_str())] -> Fill(this_pair.delta_phi, vtxZ_weight_in);
0630 
0631         if (this_pair.delta_phi > cut_GoodProtoTracklet_DeltaPhi.first && this_pair.delta_phi < cut_GoodProtoTracklet_DeltaPhi.second)
0632         {
0633             h2D_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d%s", Mbin_in, rotated_text.c_str())] -> Fill(this_pair.pair_eta_num, INTTvtxZ, vtxZ_weight_in);
0634             h2D_map[Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d%s_FineBin", Mbin_in, rotated_text.c_str())] -> Fill(this_pair.pair_eta_num, INTTvtxZ, vtxZ_weight_in);
0635         }
0636         // note : for both rotated and non-rotated
0637 
0638         // Division : ---type A---------------------------------------------------------------------------------------------------------------------------------------------------------------------    
0639         // if (ClusLadderZId->at(this_pair.inner_index) != typeA_sensorZID[0] && ClusLadderZId->at(this_pair.inner_index) != typeA_sensorZID[1]) {continue;}
0640         // if (ClusLadderZId->at(this_pair.outer_index) != typeA_sensorZID[0] && ClusLadderZId->at(this_pair.outer_index) != typeA_sensorZID[1]) {continue;}
0641         if (this_pair.inner_zid != typeA_sensorZID[0] && this_pair.inner_zid != typeA_sensorZID[1]) {continue;}
0642         if (this_pair.outer_zid != typeA_sensorZID[0] && this_pair.outer_zid != typeA_sensorZID[1]) {continue;}
0643 
0644         // note : for both rotated and non-rotated
0645         h1D_map[Form("h1D_typeA_DeltaPhi_Mbin%d_Eta%d_VtxZ%d%s", Mbin_in, eta_bin, vtxz_bin_in, rotated_text.c_str())] -> Fill(this_pair.delta_phi, vtxZ_weight_in);
0646 
0647         if (this_pair.delta_phi > cut_GoodProtoTracklet_DeltaPhi.first && this_pair.delta_phi < cut_GoodProtoTracklet_DeltaPhi.second)
0648         {
0649             h2D_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d%s", Mbin_in, rotated_text.c_str())] -> Fill(this_pair.pair_eta_num, INTTvtxZ, vtxZ_weight_in);
0650             h2D_map[Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d%s_FineBin", Mbin_in, rotated_text.c_str())] -> Fill(this_pair.pair_eta_num, INTTvtxZ, vtxZ_weight_in);
0651         }
0652         // note : for both rotated and non-rotated
0653     }
0654 
0655     if (isRotated == false){
0656         if (Pair_DeltaPhi_vec.size() != input_TrackletPair_vec.size()){
0657             std::cout<<"Something wrong: Pair_DeltaPhi_vec.size() != input_TrackletPair_vec.size()"<<std::endl;
0658             exit(1);
0659         }
0660 
0661         long long vec_size = Pair_DeltaPhi_vec.size();
0662         long long ind[Pair_DeltaPhi_vec.size()];
0663         TMath::Sort(vec_size, &Pair_DeltaPhi_vec[0], ind, false);
0664 
0665         for (int pair_i = 0; pair_i < vec_size; pair_i++){
0666 
0667             pair_str pair = input_TrackletPair_vec[ind[pair_i]];
0668             if (std::find(Used_Clus_index_vec.begin(), Used_Clus_index_vec.end(), pair.inner_index) != Used_Clus_index_vec.end()) {continue;}
0669             if (std::find(Used_Clus_index_vec.begin(), Used_Clus_index_vec.end(), pair.outer_index) != Used_Clus_index_vec.end()) {continue;}
0670 
0671             if (eID_in % 250 == 0){
0672                 std::cout<<"Event ID : "<<eID_in<<" Mbin : "<<Mbin_in<<" pair_i : "<<pair_i<<", ind[pair_i] : "<<ind[pair_i]<<", pair.delta_phi : "<<pair.delta_phi<<" pair.delta_eta : "<<pair.delta_eta<<std::endl;
0673             }
0674 
0675             if (Pair_DeltaPhi_vec[ind[pair_i]] > cut_bestPair_DeltaPhi.second) { break; }
0676 
0677             h1D_map[Form("h1D_BestPair_DeltaEta_Mbin%d_VtxZ%d",Mbin_in,vtxz_bin_in)] -> Fill(pair.delta_eta, vtxZ_weight_in);
0678             h1D_map[Form("h1D_BestPair_DeltaPhi_Mbin%d_VtxZ%d",Mbin_in,vtxz_bin_in)] -> Fill(pair.delta_phi, vtxZ_weight_in);
0679 
0680             h2D_map[Form("h2D_BestPairEtaVtxZ_Mbin%d",Mbin_in)] -> Fill(pair.pair_eta_num, INTTvtxZ, vtxZ_weight_in);
0681             h2D_map[Form("h2D_BestPairEtaVtxZ_Mbin%d_FineBin",Mbin_in)] -> Fill(pair.pair_eta_num, INTTvtxZ, vtxZ_weight_in);
0682 
0683             h1D_map[Form("h1D_BestPair_ClusPhiSize_Mbin%d_VtxZ%d",Mbin_in,vtxz_bin_in)] -> Fill(pair.inner_phi_size, vtxZ_weight_in);
0684             h1D_map[Form("h1D_BestPair_ClusPhiSize_Mbin%d_VtxZ%d",Mbin_in,vtxz_bin_in)] -> Fill(pair.outer_phi_size, vtxZ_weight_in);
0685             h1D_map[Form("h1D_BestPair_ClusAdc_Mbin%d_VtxZ%d",Mbin_in,vtxz_bin_in)] -> Fill(pair.inner_adc, vtxZ_weight_in);
0686             h1D_map[Form("h1D_BestPair_ClusAdc_Mbin%d_VtxZ%d",Mbin_in,vtxz_bin_in)] -> Fill(pair.outer_adc, vtxZ_weight_in);
0687 
0688             // note : typeA
0689             if (
0690                 // (ClusLadderZId->at(pair.inner_index) == typeA_sensorZID[0] || ClusLadderZId->at(pair.inner_index) == typeA_sensorZID[1]) &&
0691                 // (ClusLadderZId->at(pair.outer_index) == typeA_sensorZID[0] || ClusLadderZId->at(pair.outer_index) == typeA_sensorZID[1])
0692                 
0693                 (pair.inner_zid == typeA_sensorZID[0] || pair.inner_zid == typeA_sensorZID[1]) &&
0694                 (pair.outer_zid == typeA_sensorZID[0] || pair.outer_zid == typeA_sensorZID[1])
0695 
0696             ){
0697                 h1D_map[Form("h1D_typeA_BestPair_DeltaEta_Mbin%d_VtxZ%d",Mbin_in,vtxz_bin_in)] -> Fill(pair.delta_eta, vtxZ_weight_in);
0698                 h1D_map[Form("h1D_typeA_BestPair_DeltaPhi_Mbin%d_VtxZ%d",Mbin_in,vtxz_bin_in)] -> Fill(pair.delta_phi, vtxZ_weight_in);
0699 
0700                 h2D_map[Form("h2D_typeA_BestPairEtaVtxZ_Mbin%d",Mbin_in)] -> Fill(pair.pair_eta_num, INTTvtxZ, vtxZ_weight_in);
0701                 h2D_map[Form("h2D_typeA_BestPairEtaVtxZ_Mbin%d_FineBin",Mbin_in)] -> Fill(pair.pair_eta_num, INTTvtxZ, vtxZ_weight_in);
0702 
0703                 h1D_map[Form("h1D_typeA_BestPair_ClusPhiSize_Mbin%d_VtxZ%d",Mbin_in,vtxz_bin_in)] -> Fill(pair.inner_phi_size, vtxZ_weight_in);
0704                 h1D_map[Form("h1D_typeA_BestPair_ClusPhiSize_Mbin%d_VtxZ%d",Mbin_in,vtxz_bin_in)] -> Fill(pair.outer_phi_size, vtxZ_weight_in);
0705                 h1D_map[Form("h1D_typeA_BestPair_ClusAdc_Mbin%d_VtxZ%d",Mbin_in,vtxz_bin_in)] -> Fill(pair.inner_adc, vtxZ_weight_in);
0706                 h1D_map[Form("h1D_typeA_BestPair_ClusAdc_Mbin%d_VtxZ%d",Mbin_in,vtxz_bin_in)] -> Fill(pair.outer_adc, vtxZ_weight_in);
0707             }
0708 
0709 
0710             Used_Clus_index_vec.push_back(pair.inner_index);
0711             Used_Clus_index_vec.push_back(pair.outer_index);
0712         }
0713 
0714     }
0715 
0716     
0717 
0718     // h1D_PairDeltaEta_inclusive
0719     // h1D_PairDeltaPhi_inclusive
0720 
0721     // Form("h1D_DeltaPhi_Mbin%d_Eta%d_VtxZ%d", Mbin, eta_bin, vtxz_bin)
0722     // Form("h1D_typeA_DeltaPhi_Mbin%d_Eta%d_VtxZ%d", Mbin, eta_bin, vtxz_bin)
0723     // Form("h1D_DeltaPhi_Mbin%d_Eta%d_VtxZ%d_rotated", Mbin, eta_bin, vtxz_bin)
0724     // Form("h1D_typeA_DeltaPhi_Mbin%d_Eta%d_VtxZ%d_rotated", Mbin, eta_bin, vtxz_bin)
0725 
0726     // Form("h1D_BestPair_DeltaEta_Mbin%d",Mbin_in)
0727     // Form("h1D_BestPair_DeltaPhi_Mbin%d",Mbin_in)
0728     // Form("h1D_BestPair_ClusPhiSize_Mbin%d",Mbin_in)
0729     // Form("h1D_BestPair_ClusAdc_Mbin%d",Mbin_in)
0730     // Form("h1D_typeA_BestPair_DeltaEta_Mbin%d",Mbin_in)
0731     // Form("h1D_typeA_BestPair_DeltaPhi_Mbin%d",Mbin_in)
0732     // Form("h1D_typeA_BestPair_ClusPhiSize_Mbin%d",Mbin_in)
0733     // Form("h1D_typeA_BestPair_ClusAdc_Mbin%d",Mbin_in)
0734 
0735     // Form("h2D_BestPairEtaVtxZ_Mbin%d",Mbin_in)
0736     // Form("h2D_BestPairEtaVtxZ_Mbin%d_FineBin",Mbin_in)
0737     // Form("h2D_typeA_BestPairEtaVtxZ_Mbin%d",Mbin_in)
0738     // Form("h2D_typeA_BestPairEtaVtxZ_Mbin%d_FineBin",Mbin_in)
0739 
0740     // Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d", Mbin_in)
0741     // Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin", Mbin_in)
0742     // Form("h2D_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated", Mbin_in)
0743 
0744     // Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d", Mbin_in)
0745     // Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_FineBin", Mbin_in)
0746     // Form("h2D_typeA_GoodProtoTracklet_EtaVtxZ_Mbin%d_rotated", Mbin_in)
0747 
0748     
0749 }
0750 
0751 void TrackletHistogramNew::MainProcess()
0752 {
0753     if (SetRandomHits.first){PrepareUniqueClusXYZ();}
0754     
0755     if (ColMulMask && h2D_GoodColMap == nullptr){
0756         std::cout<<"The GoodColMap is not set correctly"<<std::endl;
0757         exit(1);
0758     }
0759 
0760     if (HaveGeoOffsetTag && CheckGeoOffsetMap() <= 0){
0761         std::cout<<"The geo offset map is not set correctly"<<std::endl;
0762         exit(1);
0763     }
0764 
0765     if (HaveGeoOffsetTag == false && CheckGeoOffsetMap() > 0.0001) {
0766         std::cout<<"The HaveGeoOffsetTag is false, but the GeoOffsetMap has some non-zero numbers, please check the GeoOffsetMap"<<std::endl;
0767         exit(1);
0768     }
0769 
0770     if (
0771         ColMulMask &&
0772         (
0773             h1D_GoodColMap_ZId -> GetNbinsX() != h2D_GoodColMap -> GetNbinsY() ||
0774             fabs(h1D_GoodColMap_ZId -> GetXaxis() -> GetXmin() - h2D_GoodColMap -> GetYaxis() -> GetXmin()) > 0.0000001 ||
0775             fabs(h1D_GoodColMap_ZId -> GetXaxis() -> GetXmax() - h2D_GoodColMap -> GetYaxis() -> GetXmax()) > 0.0000001
0776         )
0777 
0778     ){
0779         std::cout<<"The setting of h1D_GoodColMap_ZId is different from h2D_GoodColMap"<<std::endl;
0780         std::cout<<"h1D_GoodColMap_ZId : "<<h1D_GoodColMap_ZId -> GetNbinsX()<<" "<<h1D_GoodColMap_ZId -> GetXaxis() -> GetXmin()<<" "<<h1D_GoodColMap_ZId -> GetXaxis() -> GetXmax()<<std::endl;
0781         std::cout<<"h2D_GoodColMap : "<<h2D_GoodColMap -> GetNbinsY()<<" "<<h2D_GoodColMap -> GetYaxis() -> GetXmin()<<" "<<h2D_GoodColMap -> GetYaxis() -> GetXmax()<<std::endl;
0782         exit(1);
0783     }
0784 
0785     for (int i = 0; i < run_nEvents; i++)
0786     {
0787         tree_in -> GetEntry(i);
0788 
0789         EvtCleanUp();
0790 
0791         if (i % 10 == 0) {std::cout << "Processing event " << i<<", NClus : "<< ClusX -> size() << std::endl;}
0792 
0793         if (RandInttZ){
0794             // INTTvtxZ = rand3 -> Uniform(VtxZEdge_min,VtxZEdge_max);
0795             INTTvtxZ = rand3 -> Uniform(-10, 10);
0796             INTTvtxZError = 0.;
0797         }
0798 
0799         // =======================================================================================================================================================
0800         // note : hard cut
0801 
0802         // note : for data
0803         if (runnumber != -1 && BcoFullDiffCut && InttBcoFullDiff_next <= cut_InttBcoFullDIff_next) {continue;}
0804         if (runnumber != -1 && MBDNSg2 != 1) {continue;} // todo: assume MC no trigger
0805 
0806         // note : for MC
0807         // if (runnumber == -1 && NTruthVtx != 1) {continue;}
0808 
0809         // note : both data and MC
0810         if (is_min_bias != 1) {continue;}
0811         if (MBD_z_vtx != MBD_z_vtx) {continue;}
0812         if (MBD_centrality != MBD_centrality) {continue;}
0813         if (MBD_centrality < 0 || MBD_centrality > 100) {continue;}
0814         if (INTTvtxZ != INTTvtxZ) {continue;}
0815         if (MBD_z_vtx < cut_GlobalMBDvtxZ.first || MBD_z_vtx > cut_GlobalMBDvtxZ.second) {continue;} // todo: the hard cut 60 cm 
0816         // =======================================================================================================================================================
0817 
0818         int Mbin = h1D_centrality_bin -> Fill(MBD_centrality);
0819         Mbin = (Mbin == -1) ? -1 : Mbin - 1;
0820         if (Mbin == -1) {
0821             std::cout << "Mbin == -1, MBD_centrality = " << MBD_centrality << std::endl;
0822             continue;
0823         }
0824         // =======================================================================================================================================================
0825 
0826         // =======================================================================================================================================================
0827         // note : optional cut
0828         if (INTT_vtxZ_QA && (MBD_z_vtx - INTTvtxZ < cut_vtxZDiff.first || MBD_z_vtx - INTTvtxZ > cut_vtxZDiff.second) ) {continue;}
0829         if (INTT_vtxZ_QA && (TrapezoidalFitWidth < cut_TrapezoidalFitWidth.first || TrapezoidalFitWidth > cut_TrapezoidalFitWidth.second)){continue;}
0830         if (INTT_vtxZ_QA && (TrapezoidalFWHM < cut_TrapezoidalFWHM.first || TrapezoidalFWHM > cut_TrapezoidalFWHM.second)){continue;}
0831         if (INTT_vtxZ_QA && (INTTvtxZError < cut_INTTvtxZError.first || INTTvtxZError > cut_INTTvtxZError.second)){continue;}
0832         // =======================================================================================================================================================
0833 
0834         int InttVtxZ_bin = h1D_vtxz_template->Fill(INTTvtxZ);
0835         InttVtxZ_bin = (InttVtxZ_bin == -1) ? -1 : InttVtxZ_bin - 1;
0836 
0837         // int MBDVtxZ_bin = h1D_vtxz_template->Fill(MBD_z_vtx);
0838         // MBDVtxZ_bin = (MBDVtxZ_bin == -1) ? -1 : MBDVtxZ_bin - 1;
0839 
0840         // note : everything below is within INTT vtxZ -45 cm  ~ 45 cm
0841         if (InttVtxZ_bin == -1) {continue;}
0842 
0843 
0844         // note : for truth
0845         if (runnumber == -1){
0846             int NHadrons = 0;
0847             int NHadrons_OneEtaBin = 0;
0848 
0849             // int TruthVtxZ_bin = h1D_vtxz_template->Fill(TruthPV_trig_z);
0850             // TruthVtxZ_bin = (TruthVtxZ_bin == -1) ? -1 : TruthVtxZ_bin - 1;
0851 
0852             if (InttVtxZ_bin != -1){
0853                 for (int true_i = 0; true_i < NPrimaryG4P; true_i++){
0854                     if (PrimaryG4P_isChargeHadron->at(true_i) != 1) { continue; }
0855 
0856                     // todo : for debug
0857                     if (INTTvtxZ >= -10 && INTTvtxZ < 10){
0858                         NHadrons++;
0859                         // todo : the selected bin
0860                         if (PrimaryG4P_Eta->at(true_i) >= -1.1 && PrimaryG4P_Eta->at(true_i) < -0.9){
0861                             NHadrons_OneEtaBin++;
0862                         }
0863                     }
0864                     
0865                     h1D_map[Form("h1D_TrueEta_Mbin%d_VtxZ%d", Mbin, InttVtxZ_bin)] -> Fill(PrimaryG4P_Eta->at(true_i));
0866                     h2D_map[Form("h2D_TrueEtaVtxZ_Mbin%d", Mbin)] -> Fill(PrimaryG4P_Eta->at(true_i), INTTvtxZ);
0867                     h2D_map[Form("h2D_TrueEtaVtxZ_Mbin%d_FineBin", Mbin)] -> Fill(PrimaryG4P_Eta->at(true_i), INTTvtxZ);
0868                 }
0869 
0870                 h2D_map["h2D_TrueEvtCount_vtxZCentrality"]->Fill(INTTvtxZ, MBD_centrality);
0871 
0872                 // todo : for debug
0873                 if (INTTvtxZ >= -10 && INTTvtxZ < 10){
0874 
0875                     h1D_map["debug_h1D_NHadron_Inclusive100"] -> Fill(NHadrons);
0876                     h1D_map["debug_h1D_NHadron_OneEtaBin_Inclusive100"] -> Fill(NHadrons_OneEtaBin);
0877                 }
0878 
0879             }
0880         }
0881 
0882 
0883         // =======================================================================================================================================================
0884         
0885         if (vtxZReweight.first && runnumber != -1){
0886             std::cout<<"Should not have the vtxZ weighting from the data"<<std::endl;
0887             exit(1);
0888         }
0889 
0890         double INTTvtxZWeighting;
0891         if (vtxZReweight.first && h1D_INTT_vtxZ_reweighting != nullptr){
0892             INTTvtxZWeighting = h1D_INTT_vtxZ_reweighting -> GetBinContent(h1D_INTT_vtxZ_reweighting -> FindBin(INTTvtxZ));
0893         }
0894         else if (vtxZReweight.first && h1D_INTT_vtxZ_reweighting == nullptr){
0895             std::cout << "ApplyVtxZReWeighting is true, but h1D_INTT_vtxZ_reweighting is nullptr" << std::endl;
0896             exit(1);
0897         }
0898         else {
0899             INTTvtxZWeighting = 1.0;
0900         }
0901 
0902         // =======================================================================================================================================================
0903 
0904         // note : for reconstructed 
0905         // note : INTT vtxZ range -45 cm ~ 45 cm
0906         h1D_centrality_bin_weighted -> Fill(MBD_centrality,INTTvtxZWeighting);
0907         h1D_INTTvtxZ -> Fill(INTTvtxZ, INTTvtxZWeighting);
0908         h1D_INTTvtxZ_FineBin -> Fill(INTTvtxZ, INTTvtxZWeighting);
0909         h1D_INTTvtxZ_FineBin_NoVtxZWeight -> Fill(INTTvtxZ); // note : no vtxZ weighting
0910         h2D_INTTvtxZFineBin_CentralityBin -> Fill(INTTvtxZ, MBD_centrality, INTTvtxZWeighting);
0911         h1D_map["h1D_centrality"] -> Fill(MBD_centrality,INTTvtxZWeighting); // note : fine bins in centrality
0912 
0913         h2D_map["h2D_RecoEvtCount_vtxZCentrality"] -> Fill(INTTvtxZ, MBD_centrality, INTTvtxZWeighting);
0914         h2D_map["h2D_InttVtxZ_Centrality"] -> Fill(INTTvtxZ, MBD_centrality, INTTvtxZWeighting); // note : fine bins in centrality
0915         h1D_map[Form("h1D_centrality_InttVtxZ%d", InttVtxZ_bin)] -> Fill(MBD_centrality,INTTvtxZWeighting); // note : fine bins in centrality
0916 
0917         PrepareClusterVec();
0918 
0919         for (ClusHistogram::clu_info this_clu : evt_sPH_inner_nocolumn_vec){
0920             h1D_Inner_ClusEta_INTTz -> Fill(this_clu.eta_INTTz, INTTvtxZWeighting);
0921         }
0922         for (ClusHistogram::clu_info this_clu : evt_sPH_outer_nocolumn_vec){
0923             h1D_Outer_ClusEta_INTTz -> Fill(this_clu.eta_INTTz, INTTvtxZWeighting);
0924         }
0925         
0926 
0927         GetTrackletPair(evt_TrackletPair_vec, false);
0928         FillPairs(evt_TrackletPair_vec, false, Mbin, InttVtxZ_bin, INTTvtxZWeighting, i);
0929 
0930         GetTrackletPair(evt_TrackletPairRotate_vec, true);  
0931         FillPairs(evt_TrackletPairRotate_vec, true, Mbin, InttVtxZ_bin, INTTvtxZWeighting, i);
0932 
0933     }
0934 }
0935 
0936 void TrackletHistogramNew::EndRun()
0937 {
0938     file_out -> cd();
0939 
0940     tree_out_par -> Fill();
0941     tree_out_par -> Write();
0942 
0943     h1D_INTTvtxZ -> Write();
0944     h1D_INTTvtxZ_FineBin -> Write();
0945     h1D_INTTvtxZ_FineBin_NoVtxZWeight -> Write();
0946     h2D_INTTvtxZFineBin_CentralityBin -> Write();
0947 
0948     h1D_centrality_bin -> Write();
0949     h1D_centrality_bin_weighted -> Write();
0950     
0951     h1D_eta_template -> Write();
0952     h1D_vtxz_template -> Write();
0953 
0954     h1D_PairDeltaEta_inclusive -> Write();
0955     h1D_PairDeltaPhi_inclusive -> Write();
0956 
0957     h1D_Inner_ClusEta_INTTz -> Write();
0958     h1D_Outer_ClusEta_INTTz -> Write();
0959 
0960     if (h1D_RandClusZ_ref != nullptr){
0961         h1D_RandClusZ_ref -> Write();
0962     }
0963 
0964     if (h2D_RandClusXY_ref != nullptr){
0965         h2D_RandClusXY_ref -> Write();
0966     }
0967 
0968     for (auto &pair : h2D_map){
0969         pair.second -> Write();
0970     }
0971 
0972     for (auto &pair : h1D_map){
0973         pair.second -> Write();
0974     }
0975 
0976     file_out -> Close();
0977 }
0978 
0979 std::vector<ClusHistogram::clu_info> TrackletHistogramNew::GetRotatedClusterVec(std::vector<ClusHistogram::clu_info> input_cluster_vec)
0980 {
0981     std::vector<ClusHistogram::clu_info> output_cluster_vec; output_cluster_vec.clear();
0982 
0983     for (ClusHistogram::clu_info this_clu : input_cluster_vec)
0984     {
0985         std::pair<double,double> rotated_xy = rotatePoint(this_clu.x, this_clu.y);
0986      
0987         ClusHistogram::clu_info rotated_clu = this_clu;
0988         rotated_clu.x = rotated_xy.first;
0989         rotated_clu.y = rotated_xy.second;
0990 
0991         rotated_clu.eta_INTTz = ClusHistogram::get_clu_eta({vertexXYIncm.first, vertexXYIncm.second, INTTvtxZ}, {rotated_clu.x, rotated_clu.y, rotated_clu.z});
0992 
0993         output_cluster_vec.push_back(rotated_clu);
0994     }
0995 
0996     return output_cluster_vec;
0997 }
0998 
0999 std::pair<double,double> TrackletHistogramNew::rotatePoint(double x, double y)
1000 {
1001     // Convert the rotation angle from degrees to radians
1002     double rotation = rotate_phi_angle; // todo rotation is here
1003     double angleRad = rotation * M_PI / 180.0;
1004 
1005     // Perform the rotation
1006     double xOut = x * cos(angleRad) - y * sin(angleRad);
1007     double yOut = x * sin(angleRad) + y * cos(angleRad);
1008 
1009     xOut = (fabs(xOut) < 0.0000001) ? 0 : xOut;
1010     yOut = (fabs(yOut) < 0.0000001) ? 0 : yOut;
1011 
1012     // cout<<"Post rotation: "<<xOut<<" "<<yOut<<endl;
1013 
1014     return {xOut,yOut};
1015 }
1016 
1017 std::pair<double,double> TrackletHistogramNew::Get_possible_zvtx(double rvtx, std::vector<double> p0, std::vector<double> p1) // note : inner p0, outer p1, vector {r,z, zid}, -> {y,x}
1018 {
1019     std::vector<double> p0_z_edge = { ( p0[2] == typeA_sensorZID[0] || p0[2] == typeA_sensorZID[1] ) ? p0[1] - typeA_sensor_half_length_incm : p0[1] - typeB_sensor_half_length_incm, ( p0[2] == typeA_sensorZID[0] || p0[2] == typeA_sensorZID[1] ) ? p0[1] + typeA_sensor_half_length_incm : p0[1] + typeB_sensor_half_length_incm}; // note : vector {left edge, right edge}
1020     std::vector<double> p1_z_edge = { ( p1[2] == typeA_sensorZID[0] || p1[2] == typeA_sensorZID[1] ) ? p1[1] - typeA_sensor_half_length_incm : p1[1] - typeB_sensor_half_length_incm, ( p1[2] == typeA_sensorZID[0] || p1[2] == typeA_sensorZID[1] ) ? p1[1] + typeA_sensor_half_length_incm : p1[1] + typeB_sensor_half_length_incm}; // note : vector {left edge, right edge}
1021 
1022     double edge_first  = Get_extrapolation(rvtx,p0_z_edge[0],p0[0],p1_z_edge[1],p1[0]);
1023     double edge_second = Get_extrapolation(rvtx,p0_z_edge[1],p0[0],p1_z_edge[0],p1[0]);
1024 
1025     double mid_point = (edge_first + edge_second) / 2.;
1026     double possible_width = fabs(edge_first - edge_second) / 2.;
1027 
1028     return {mid_point, possible_width}; // note : first : mid point, second : width
1029 
1030 }
1031 
1032 // note : angle_1 = inner clu phi
1033 // note: angle_2 = outer clu phi
1034 double TrackletHistogramNew::get_delta_phi(double angle_1, double angle_2)
1035 {
1036     std::vector<double> vec_abs = {std::fabs(angle_1 - angle_2), std::fabs(angle_1 - angle_2 + 360), std::fabs(angle_1 - angle_2 - 360)};
1037     std::vector<double> vec = {(angle_1 - angle_2), (angle_1 - angle_2 + 360), (angle_1 - angle_2 - 360)};
1038     return vec[std::distance(vec_abs.begin(), std::min_element(vec_abs.begin(),vec_abs.end()))];
1039 }
1040 
1041 double TrackletHistogramNew::get_radius(double x, double y)
1042 {
1043     return sqrt(pow(x,2)+pow(y,2));
1044 }
1045 
1046 // note : pair<reduced chi2, eta of the track>
1047 // note : vector : {r,z}
1048 // note : p0 : vertex point {r,z,zerror}
1049 // note : p1 : inner layer
1050 // note : p2 : outer layer
1051 std::pair<double,double> TrackletHistogramNew::Get_eta(std::vector<double>p0, std::vector<double>p1, std::vector<double>p2)
1052 {
1053     // if (track_gr -> GetN() != 0) {cout<<"In TrackletHistogramNew class, track_gr is not empty, track_gr size : "<<track_gr -> GetN()<<endl; exit(0);}
1054     
1055     if (track_gr -> GetN() != 0) {track_gr -> Set(0);}
1056     
1057     std::vector<double> r_vec  = {p0[0], p1[0], p2[0]}; 
1058     std::vector<double> z_vec  = {p0[1], p1[1], p2[1]}; 
1059     std::vector<double> re_vec = {0, 0, 0}; 
1060     std::vector<double> ze_vec = {
1061         p0[2],
1062         (p1[2] == typeA_sensorZID[0] || p1[2] == typeA_sensorZID[1]) ? typeA_sensor_half_length_incm : typeB_sensor_half_length_incm,
1063         (p2[2] == typeA_sensorZID[0] || p2[2] == typeA_sensorZID[1]) ? typeA_sensor_half_length_incm : typeB_sensor_half_length_incm
1064     };
1065 
1066     // std::vector<double> ze_vec = {p0[2], ( fabs( p1[1] ) < sensor_type_segment ) ? sensor_length_typeA : sensor_length_typeB, ( fabs( p2[1] ) < sensor_type_segment ) ? sensor_length_typeA : sensor_length_typeB}; 
1067 
1068     // note : swap the r and z, to have easier fitting 
1069     // note : in principle, Z should be in the x axis, R should be in the Y axis
1070     for (int i = 0; i < 3; i++)
1071     {
1072         track_gr -> SetPoint(i,r_vec[i],z_vec[i]);
1073         track_gr -> SetPointError(i,re_vec[i],ze_vec[i]);
1074 
1075         // cout<<"In TrackletHistogramNew class, point : "<<i<<" r : "<<r_vec[i]<<" z : "<<z_vec[i]<<" re : "<<re_vec[i]<<" ze : "<<ze_vec[i]<<endl;
1076     }    
1077     
1078     double vertical_line = ( fabs( grEY_stddev(track_gr) ) < 0.00001 ) ? 1 : 0;
1079     
1080     if (vertical_line) {return {0,0};}
1081     else 
1082     {
1083         fit_rz -> SetParameters(0,0);
1084 
1085         track_gr -> Fit(fit_rz,"NQ");
1086 
1087         std::pair<double,double> ax_b = mirrorPolynomial(fit_rz -> GetParameter(1),fit_rz -> GetParameter(0));
1088 
1089         return  {(fit_rz -> GetChisquare() / double(fit_rz -> GetNDF())), -1 * TMath::Log( fabs( tan( atan2(ax_b.first, (ax_b.first > 0) ? 1. : -1. ) / 2 ) ) )};
1090     }
1091 
1092 }
1093 
1094 double TrackletHistogramNew::grEY_stddev(TGraphErrors * input_grr)
1095 {
1096     std::vector<double> input_vector; input_vector.clear();
1097     for (int i = 0; i < input_grr -> GetN(); i++) {input_vector.push_back( input_grr -> GetPointY(i) );}
1098 
1099     double average = std::accumulate( input_vector.begin(), input_vector.end(), 0.0 ) / double(input_vector.size());
1100 
1101     double sum_subt = 0;
1102     for (int ele : input_vector) {sum_subt+=pow((ele-average),2);}
1103     
1104     return sqrt(sum_subt/(input_vector.size()-1));
1105 }   
1106 
1107 std::pair<double, double> TrackletHistogramNew::mirrorPolynomial(double a, double b) {
1108     // Interchange 'x' and 'y'
1109     double mirroredA = 1.0 / a;
1110     double mirroredB = -b / a;
1111 
1112     return {mirroredA, mirroredB};
1113 }
1114 
1115 double TrackletHistogramNew::Get_extrapolation(double given_y, double p0x, double p0y, double p1x, double p1y) // note : x : z, y : r
1116 {
1117     if ( fabs(p0x - p1x) < 0.000001 ){ // note : the line is vertical (if z is along the x axis)
1118         return p0x;
1119     }
1120     else {
1121         double slope = (p1y - p0y) / (p1x - p0x);
1122         double yIntercept = p0y - slope * p0x;
1123         double xCoordinate = (given_y - yIntercept) / slope;
1124         return xCoordinate;
1125     }
1126 }
[ 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!