Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "ClusHistogram.h"
0002 
0003 ClusHistogram::ClusHistogram(
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, // note : {adc, phi size}
0018         bool HaveGeoOffsetTag_in,
0019         std::pair<bool, int> SetRandomHits_in,
0020         bool RandInttZ_in,
0021         bool ColMulMask_in,
0022         int c_type_in // note : constructor type
0023 ):
0024     process_id(process_id_in),
0025     runnumber(runnumber_in),
0026     run_nEvents(run_nEvents_in),
0027     input_directory(input_directory_in),
0028     input_file_name(input_file_name_in),
0029     output_directory(output_directory_in),
0030 
0031     output_file_name_suffix(output_file_name_suffix_in),
0032     vertexXYIncm(vertexXYIncm_in),
0033 
0034     vtxZReweight(vtxZReweight_in),
0035     BcoFullDiffCut(BcoFullDiffCut_in),
0036     INTT_vtxZ_QA(INTT_vtxZ_QA_in),
0037     isClusQA(isClusQA_in),
0038     HaveGeoOffsetTag(HaveGeoOffsetTag_in),
0039     SetRandomHits(SetRandomHits_in),
0040     RandInttZ(RandInttZ_in),
0041     ColMulMask(ColMulMask_in),
0042 
0043     c_type(c_type_in) // note : the constructor switch
0044 {
0045     PrepareInputRootFile();
0046 
0047     h1D_INTT_vtxZ_reweighting = (vtxZReweight.first) ? (TH1D*)vtxZReweight.second->Clone() : nullptr;
0048 
0049     nCentrality_bin = centrality_edges.size() - 1;
0050 
0051     run_nEvents = (run_nEvents == -1) ? tree_in->GetEntries() : run_nEvents;
0052     run_nEvents = (run_nEvents > tree_in->GetEntries()) ? tree_in->GetEntries() : run_nEvents;
0053 
0054     if (HaveGeoOffsetTag == false){
0055         for (int layer_i = B0L0_index; layer_i <= B1L1_index; layer_i++){
0056             double N_layer_ladder = (layer_i == B0L0_index || layer_i == B0L1_index) ? nLadder_inner : nLadder_outer;
0057             for (int phi_i = 0; phi_i < N_layer_ladder; phi_i++){
0058                 geo_offset_map[Form("%i_%i", layer_i, phi_i)] = {0, 0, 0};
0059             }
0060         }
0061     }
0062 
0063     if (RandInttZ && INTT_vtxZ_QA){
0064         std::cout << "RandInttZ and INTT_vtxZ_QA cannot be true at the same time" << std::endl;
0065         exit(1);
0066     }
0067 
0068     if (SetRandomHits.first || RandInttZ){
0069         rand3 = new TRandom3(0);
0070     }
0071 
0072     if (c_type == 0){
0073         PrepareOutPutFileName();
0074         PrepareOutPutRootFile();
0075         PrepareHistograms();
0076     }
0077 
0078 }
0079 
0080 std::map<std::string, int> ClusHistogram::GetInputTreeBranchesMap(TTree * m_tree_in)
0081 {
0082     std::map<std::string, int> branch_map;
0083     TObjArray * branch_list = m_tree_in -> GetListOfBranches();
0084     for (int i = 0; i < branch_list -> GetEntries(); i++)
0085     {
0086         TBranch * branch = dynamic_cast<TBranch*>(branch_list->At(i));
0087         branch_map[branch -> GetName()] = 1;
0088     }
0089     return branch_map;
0090 }
0091 
0092 void ClusHistogram::PrepareInputRootFile()
0093 {
0094     file_in = TFile::Open(Form("%s/%s", input_directory.c_str(), input_file_name.c_str()));
0095     if (!file_in || file_in -> IsZombie() || file_in == nullptr) {
0096         std::cout << "Error: cannot open file: " << input_file_name << std::endl;
0097         exit(1);
0098     }
0099 
0100     tree_in = (TTree*)file_in -> Get("EventTree");
0101     
0102     std::map<std::string, int> branch_map = GetInputTreeBranchesMap(tree_in);
0103     if(branch_map.find("event") != branch_map.end()){tree_in -> SetBranchStatus("event",0);}
0104 
0105     // if(branch_map.find("ClusEta_INTTz") != branch_map.end()){tree_in -> SetBranchStatus("ClusEta_INTTz",0);}
0106     if(branch_map.find("ClusEta_MBDz") != branch_map.end()){tree_in -> SetBranchStatus("ClusEta_MBDz",0);}
0107     if(branch_map.find("ClusPhi_AvgPV") != branch_map.end()){tree_in -> SetBranchStatus("ClusPhi_AvgPV",0);}
0108     if(branch_map.find("ClusEta_TrueXYZ") != branch_map.end()){tree_in -> SetBranchStatus("ClusEta_TrueXYZ",0);}
0109     if(branch_map.find("ClusPhi_TrueXY") != branch_map.end()){tree_in -> SetBranchStatus("ClusPhi_TrueXY",0);}
0110 
0111     
0112     if(branch_map.find("INTT_BCO") != branch_map.end()){tree_in -> SetBranchStatus("INTT_BCO",0);}
0113     if(branch_map.find("ClusTrkrHitSetKey") != branch_map.end()){tree_in -> SetBranchStatus("ClusTrkrHitSetKey",0);}
0114     if(branch_map.find("ClusTimeBucketId") != branch_map.end()){tree_in -> SetBranchStatus("ClusTimeBucketId",0);}
0115     
0116     if(branch_map.find("GL1Packet_BCO") != branch_map.end()){tree_in -> SetBranchStatus("GL1Packet_BCO",0);}
0117     if(branch_map.find("ncoll") != branch_map.end()){tree_in -> SetBranchStatus("ncoll",0);}
0118     if(branch_map.find("npart") != branch_map.end()){tree_in -> SetBranchStatus("npart",0);}
0119     if(branch_map.find("centrality_bimp") != branch_map.end()){tree_in -> SetBranchStatus("centrality_bimp",0);}
0120     if(branch_map.find("centrality_impactparam") != branch_map.end()){tree_in -> SetBranchStatus("centrality_impactparam",0);}
0121     if(branch_map.find("clk") != branch_map.end()){tree_in -> SetBranchStatus("clk",0);}
0122     if(branch_map.find("femclk") != branch_map.end()){tree_in -> SetBranchStatus("femclk",0);}
0123     if(branch_map.find("is_min_bias_wozdc") != branch_map.end()){tree_in -> SetBranchStatus("is_min_bias_wozdc",0);}
0124     if(branch_map.find("MBD_south_npmt") != branch_map.end()){tree_in -> SetBranchStatus("MBD_south_npmt",0);}
0125     if(branch_map.find("MBD_north_npmt") != branch_map.end()){tree_in -> SetBranchStatus("MBD_north_npmt",0);}
0126     if(branch_map.find("MBD_nhitsoverths_south") != branch_map.end()){tree_in -> SetBranchStatus("MBD_nhitsoverths_south",0);}
0127     if(branch_map.find("MBD_nhitsoverths_north") != branch_map.end()){tree_in -> SetBranchStatus("MBD_nhitsoverths_north",0);}
0128     
0129     if(branch_map.find("TrackletPair") != branch_map.end()) {tree_in -> SetBranchStatus("TrackletPair", 0);}
0130     if(branch_map.find("TrackletPairRotate") != branch_map.end()) {tree_in -> SetBranchStatus("TrackletPairRotate", 0);}
0131 
0132 
0133     // note: for reading 
0134     ClusX = 0;
0135     ClusY = 0;
0136     ClusZ = 0;
0137     ClusLayer = 0;
0138     ClusLadderZId = 0;
0139     ClusLadderPhiId = 0;
0140     ClusAdc = 0;
0141     ClusPhiSize = 0;
0142     ClusEta_INTTz = 0;
0143     
0144     // evt_TrackletPair_vec = 0;
0145     // evt_TrackletPairRotate_vec = 0;
0146 
0147     PrimaryG4P_Pt = 0;
0148     PrimaryG4P_Eta = 0;
0149     PrimaryG4P_Phi = 0;
0150     PrimaryG4P_E = 0;
0151     PrimaryG4P_PID = 0;
0152     PrimaryG4P_isChargeHadron = 0;
0153 
0154     tree_in -> SetBranchAddress("MBD_z_vtx", &MBD_z_vtx);
0155     tree_in -> SetBranchAddress("is_min_bias", &is_min_bias);
0156     tree_in -> SetBranchAddress("MBD_centrality", &MBD_centrality);
0157     tree_in -> SetBranchAddress("MBD_south_charge_sum", &MBD_south_charge_sum);
0158     tree_in -> SetBranchAddress("MBD_north_charge_sum", &MBD_north_charge_sum);
0159     tree_in -> SetBranchAddress("MBD_charge_sum", &MBD_charge_sum);
0160     tree_in -> SetBranchAddress("MBD_charge_asymm", &MBD_charge_asymm);
0161     if (runnumber != -1) {
0162         if(branch_map.find("InttBcoFullDiff_next") != branch_map.end()){
0163             tree_in -> SetBranchAddress("InttBcoFullDiff_next", &InttBcoFullDiff_next);
0164         }
0165     }
0166 
0167     // note : trigger tag
0168     if (runnumber != -1){ // note : for data
0169         if(branch_map.find("MBDNSg2") != branch_map.end()) {tree_in -> SetBranchAddress("MBDNSg2", &MBDNSg2);}
0170         if(branch_map.find("MBDNSg2_vtxZ10cm") != branch_map.end()) {tree_in -> SetBranchAddress("MBDNSg2_vtxZ10cm", &MBDNSg2_vtxZ10cm);}
0171         if(branch_map.find("MBDNSg2_vtxZ30cm") != branch_map.end()) {tree_in -> SetBranchAddress("MBDNSg2_vtxZ30cm", &MBDNSg2_vtxZ30cm);}
0172         if(branch_map.find("MBDNSg2_vtxZ60cm") != branch_map.end()) {tree_in -> SetBranchAddress("MBDNSg2_vtxZ60cm", &MBDNSg2_vtxZ60cm);}
0173     }
0174 
0175     tree_in -> SetBranchAddress("ClusX", &ClusX);
0176     tree_in -> SetBranchAddress("ClusY", &ClusY);
0177     tree_in -> SetBranchAddress("ClusZ", &ClusZ);
0178     tree_in -> SetBranchAddress("ClusLayer", &ClusLayer);
0179     tree_in -> SetBranchAddress("ClusLadderZId", &ClusLadderZId);
0180     tree_in -> SetBranchAddress("ClusLadderPhiId", &ClusLadderPhiId);
0181     tree_in -> SetBranchAddress("ClusAdc", &ClusAdc);
0182     tree_in -> SetBranchAddress("ClusPhiSize", &ClusPhiSize);
0183     if(branch_map.find("ClusEta_INTTz") != branch_map.end()) {tree_in -> SetBranchAddress("ClusEta_INTTz", &ClusEta_INTTz);}
0184 
0185     // note : INTT vertex Z
0186     if(branch_map.find("INTTvtxZ") != branch_map.end()) {tree_in -> SetBranchAddress("INTTvtxZ", &INTTvtxZ);}
0187     if(branch_map.find("INTTvtxZError") != branch_map.end()) {tree_in -> SetBranchAddress("INTTvtxZError", &INTTvtxZError);}
0188     if(branch_map.find("NgroupTrapezoidal") != branch_map.end()) {tree_in -> SetBranchAddress("NgroupTrapezoidal", &NgroupTrapezoidal);}
0189     if(branch_map.find("NgroupCoarse") != branch_map.end()) {tree_in -> SetBranchAddress("NgroupCoarse", &NgroupCoarse);}
0190     if(branch_map.find("TrapezoidalFitWidth") != branch_map.end()) {tree_in -> SetBranchAddress("TrapezoidalFitWidth", &TrapezoidalFitWidth);}
0191     if(branch_map.find("TrapezoidalFWHM") != branch_map.end()) {tree_in -> SetBranchAddress("TrapezoidalFWHM", &TrapezoidalFWHM);}
0192 
0193     // note : the tracklet pair
0194     // tree_in -> SetBranchAddress("TrackletPair", &evt_TrackletPair_vec);
0195     // tree_in -> SetBranchAddress("TrackletPairRotate", &evt_TrackletPairRotate_vec);
0196 
0197     // note : MC
0198     if (runnumber == -1){
0199         tree_in -> SetBranchAddress("TruthPV_trig_x", &TruthPV_trig_x);
0200         tree_in -> SetBranchAddress("TruthPV_trig_y", &TruthPV_trig_y);
0201         tree_in -> SetBranchAddress("TruthPV_trig_z", &TruthPV_trig_z);
0202         tree_in -> SetBranchAddress("NTruthVtx", &NTruthVtx);
0203         tree_in -> SetBranchAddress("NPrimaryG4P", &NPrimaryG4P);
0204         tree_in -> SetBranchAddress("PrimaryG4P_Pt", &PrimaryG4P_Pt);
0205         tree_in -> SetBranchAddress("PrimaryG4P_Eta", &PrimaryG4P_Eta);
0206         tree_in -> SetBranchAddress("PrimaryG4P_Phi", &PrimaryG4P_Phi);
0207         tree_in -> SetBranchAddress("PrimaryG4P_E", &PrimaryG4P_E);
0208         tree_in -> SetBranchAddress("PrimaryG4P_PID", &PrimaryG4P_PID);
0209         tree_in -> SetBranchAddress("PrimaryG4P_isChargeHadron", &PrimaryG4P_isChargeHadron);
0210     }
0211 }
0212 
0213 void ClusHistogram::PrepareOutPutFileName()
0214 {
0215     std::string job_index = std::to_string( process_id );
0216     int job_index_len = 5;
0217     job_index.insert(0, job_index_len - job_index.size(), '0');
0218 
0219     std::string runnumber_str = std::to_string( runnumber );
0220     if (runnumber != -1){
0221         int runnumber_str_len = 8;
0222         runnumber_str.insert(0, runnumber_str_len - runnumber_str.size(), '0');
0223     }
0224 
0225     if (output_file_name_suffix.size() > 0 && output_file_name_suffix[0] != '_') {
0226         output_file_name_suffix = "_" + output_file_name_suffix;
0227     }
0228 
0229     output_filename = "ClusHist";
0230     output_filename = (runnumber != -1) ? "Data_" + output_filename : "MC_" + output_filename;
0231     output_filename += (vtxZReweight.first) ? "_vtxZReweight" : "";
0232     output_filename += (BcoFullDiffCut && runnumber != -1) ? "_BcoFullDiffCut" : "";
0233     output_filename += (INTT_vtxZ_QA) ? "_VtxZQA" : "";
0234     output_filename += (isClusQA.first) ? Form("_ClusQAAdc%.0fPhiSize%.0f",isClusQA.second.first,isClusQA.second.second) : "";
0235     output_filename += (ColMulMask) ? "_ColMulMask" : "";
0236     output_filename += (HaveGeoOffsetTag) ? "_GeoOffset" : "";
0237     output_filename += (SetRandomHits.first) ? Form("_Rand%dHits",SetRandomHits.second) : "";
0238     output_filename += (RandInttZ) ? "_RandInttZ" : "";
0239     output_filename += output_file_name_suffix;
0240     output_filename += (runnumber != -1) ? Form("_%s_%s.root",runnumber_str.c_str(),job_index.c_str()) : Form("_%s.root",job_index.c_str());
0241 }
0242 
0243 
0244 void ClusHistogram::PrepareOutPutRootFile()
0245 {
0246     file_out = new TFile(Form("%s/%s", output_directory.c_str(), output_filename.c_str()), "RECREATE");
0247     tree_out_par = new TTree("tree_par", "used parameters");
0248 
0249     tree_out_par -> Branch("process_id", &process_id);
0250     tree_out_par -> Branch("runnumber", &runnumber);
0251     tree_out_par -> Branch("vtxZReweight", &vtxZReweight.first);
0252     tree_out_par -> Branch("BcoFullDiffCut", &BcoFullDiffCut);
0253     tree_out_par -> Branch("INTT_vtxZ_QA", &INTT_vtxZ_QA);
0254     tree_out_par -> Branch("isClusQA", &isClusQA.first);
0255     tree_out_par -> Branch("isClusQA_adc", &isClusQA.second.first);
0256     tree_out_par -> Branch("isClusQA_phiSize", &isClusQA.second.second);
0257 
0258     tree_out_par -> Branch("centrality_edges", &centrality_edges);
0259     tree_out_par -> Branch("nCentrality_bin", &nCentrality_bin);
0260 
0261     tree_out_par -> Branch("CentralityFineEdge_min", &CentralityFineEdge_min);
0262     tree_out_par -> Branch("CentralityFineEdge_max", &CentralityFineEdge_max);
0263     tree_out_par -> Branch("nCentralityFineBin", &nCentralityFineBin);
0264 
0265     tree_out_par -> Branch("EtaEdge_min", &EtaEdge_min);
0266     tree_out_par -> Branch("EtaEdge_max", &EtaEdge_max);
0267     tree_out_par -> Branch("nEtaBin", &nEtaBin);
0268 
0269     tree_out_par -> Branch("VtxZEdge_min", &VtxZEdge_min);
0270     tree_out_par -> Branch("VtxZEdge_max", &VtxZEdge_max);
0271     tree_out_par -> Branch("nVtxZBin", &nVtxZBin);
0272 
0273     tree_out_par -> Branch("typeA_sensorZID", &typeA_sensorZID);
0274     tree_out_par -> Branch("cut_GlobalMBDvtxZ", &cut_GlobalMBDvtxZ);   
0275     tree_out_par -> Branch("cut_vtxZDiff", &cut_vtxZDiff);
0276     tree_out_par -> Branch("cut_TrapezoidalFitWidth", &cut_TrapezoidalFitWidth);
0277     tree_out_par -> Branch("cut_TrapezoidalFWHM", &cut_TrapezoidalFWHM);
0278     tree_out_par -> Branch("cut_INTTvtxZError", &cut_INTTvtxZError);
0279 
0280     tree_out_par -> Branch("cut_InttBcoFullDIff_next", &cut_InttBcoFullDIff_next);
0281 }
0282 
0283 void ClusHistogram::PrepareHistograms()
0284 {
0285     h1D_INTTvtxZ = new TH1D("h1D_INTTvtxZ", "h1D_INTTvtxZ;INTT vtxZ [cm];Entries", 60, -60, 60);
0286 
0287     h1D_centrality_bin = new TH1D("h1D_centrality_bin","h1D_centrality_bin;Centrality [%];Entries",nCentrality_bin,&centrality_edges[0]); // note : the 0-5%
0288     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%
0289 
0290     h1D_vtxz_template = new TH1D("h1D_vtxz_template", "h1D_vtxz_template", nVtxZBin, VtxZEdge_min, VtxZEdge_max); // note : coarse
0291 
0292     h1D_GoodColMap_ZId = new TH1D("h1D_GoodColMap_ZId","h1D_GoodColMap_ZId;ClusZ [cm];Entries",nZbin, Zmin, Zmax);
0293 
0294     h1D_centrality_map.clear();
0295     h1D_centrality_map.insert( std::make_pair(
0296             "h1D_centrality",
0297             new TH1D("h1D_centrality", "h1D_centrality;Centrality [%];Entries", nCentralityFineBin, CentralityFineEdge_min, CentralityFineEdge_max)
0298         )
0299     );
0300 
0301     for (int VtxZ_bin = 0; VtxZ_bin < nVtxZBin; VtxZ_bin++)
0302     {
0303         h1D_centrality_map.insert( std::make_pair(
0304                 Form("h1D_centrality_InttVtxZ%d", VtxZ_bin),
0305                 new TH1D(Form("h1D_centrality_InttVtxZ%d", VtxZ_bin), Form("h1D_centrality_InttVtxZ%d;Centrality [%];Entries",VtxZ_bin), nCentralityFineBin, CentralityFineEdge_min, CentralityFineEdge_max)
0306             )
0307         );
0308     }
0309 
0310     h1D_NClusEta_map.clear(); // note : {inner, outer, typeA, all} x {Mbin x vtxZ bin}
0311     for (int VtxZ_bin = 0; VtxZ_bin < nVtxZBin; VtxZ_bin++){
0312         for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
0313 
0314             h1D_NClusEta_map.insert( // note : all, inner
0315                 std::make_pair(
0316                     Form("h1D_inner_NClusEta_Mbin%d_VtxZ%d", Mbin, VtxZ_bin),
0317                     new TH1D(Form("h1D_inner_NClusEta_Mbin%d_VtxZ%d", Mbin, VtxZ_bin), Form("h1D_inner_NClusEta_Mbin%d_VtxZ%d; Clus #eta (inner);Entries", Mbin, VtxZ_bin), nEtaBin, EtaEdge_min, EtaEdge_max)
0318                 )
0319             );
0320 
0321             h1D_NClusEta_map.insert( // note : all, outer
0322                 std::make_pair(
0323                     Form("h1D_outer_NClusEta_Mbin%d_VtxZ%d", Mbin, VtxZ_bin),
0324                     new TH1D(Form("h1D_outer_NClusEta_Mbin%d_VtxZ%d", Mbin, VtxZ_bin), Form("h1D_outer_NClusEta_Mbin%d_VtxZ%d; Clus #eta (outer);Entries", Mbin, VtxZ_bin), nEtaBin, EtaEdge_min, EtaEdge_max)
0325                 )
0326             );
0327 
0328             h1D_NClusEta_map.insert( // note : typeA, inner
0329                 std::make_pair(
0330                     Form("h1D_typeA_inner_NClusEta_Mbin%d_VtxZ%d", Mbin, VtxZ_bin),
0331                     new TH1D(Form("h1D_typeA_inner_NClusEta_Mbin%d_VtxZ%d", Mbin, VtxZ_bin), Form("h1D_typeA_inner_NClusEta_Mbin%d_VtxZ%d; Clus #eta (inner, type A);Entries", Mbin, VtxZ_bin), nEtaBin, EtaEdge_min, EtaEdge_max)
0332                 )
0333             );
0334 
0335             h1D_NClusEta_map.insert( // note : typeA, outer
0336                 std::make_pair(
0337                     Form("h1D_typeA_outer_NClusEta_Mbin%d_VtxZ%d", Mbin, VtxZ_bin),
0338                     new TH1D(Form("h1D_typeA_outer_NClusEta_Mbin%d_VtxZ%d", Mbin, VtxZ_bin), Form("h1D_typeA_outer_NClusEta_Mbin%d_VtxZ%d; Clus #eta (outer, type A);Entries", Mbin, VtxZ_bin), nEtaBin, EtaEdge_min, EtaEdge_max)
0339                 )
0340             );
0341         }
0342     }
0343 
0344     h2D_NClusEtaVtxZ_map.clear(); // note : {inner, outer, typeA, all} x {Mbin}
0345     for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++){
0346 
0347         h2D_NClusEtaVtxZ_map.insert( // note : all, inner
0348             std::make_pair(
0349                 Form("h2D_inner_NClusEta_Mbin%d", Mbin),
0350                 new TH2D(Form("h2D_inner_NClusEta_Mbin%d", Mbin), Form("h2D_inner_NClusEta_Mbin%d; Clus #eta (inner);INTT vtxZ [cm]", Mbin), nEtaBin, EtaEdge_min, EtaEdge_max, nVtxZBin, VtxZEdge_min, VtxZEdge_max)
0351             )
0352         );
0353 
0354         h2D_NClusEtaVtxZ_map.insert( // note : all, outer
0355             std::make_pair(
0356                 Form("h2D_outer_NClusEta_Mbin%d", Mbin),
0357                 new TH2D(Form("h2D_outer_NClusEta_Mbin%d", Mbin), Form("h2D_outer_NClusEta_Mbin%d; Clus #eta (outer);INTT vtxZ [cm]", Mbin), nEtaBin, EtaEdge_min, EtaEdge_max, nVtxZBin, VtxZEdge_min, VtxZEdge_max)
0358             )
0359         );
0360 
0361         h2D_NClusEtaVtxZ_map.insert( // note : typeA, inner
0362             std::make_pair(
0363                 Form("h2D_typeA_inner_NClusEta_Mbin%d", Mbin),
0364                 new TH2D(Form("h2D_typeA_inner_NClusEta_Mbin%d", Mbin), Form("h2D_typeA_inner_NClusEta_Mbin%d; Clus #eta (inner, type A);INTT vtxZ [cm]", Mbin), nEtaBin, EtaEdge_min, EtaEdge_max, nVtxZBin, VtxZEdge_min, VtxZEdge_max)
0365             )
0366         );
0367 
0368         h2D_NClusEtaVtxZ_map.insert( // note : typeA, outer
0369             std::make_pair(
0370                 Form("h2D_typeA_outer_NClusEta_Mbin%d", Mbin),
0371                 new TH2D(Form("h2D_typeA_outer_NClusEta_Mbin%d", Mbin), Form("h2D_typeA_outer_NClusEta_Mbin%d; Clus #eta (outer, type A);INTT vtxZ [cm]", Mbin), nEtaBin, EtaEdge_min, EtaEdge_max, nVtxZBin, VtxZEdge_min, VtxZEdge_max)
0372             )
0373         );
0374 
0375         // Division: ---Fine Bin-----------------------------------------------------------------------------------------------------------------------------------------
0376         h2D_NClusEtaVtxZ_map.insert( // note : all, inner
0377             std::make_pair(
0378                 Form("h2D_inner_NClusEta_Mbin%d_FineBin", Mbin),
0379                 new TH2D(Form("h2D_inner_NClusEta_Mbin%d_FineBin", Mbin), Form("h2D_inner_NClusEta_Mbin%d_FineBin; Clus #eta (inner);INTT vtxZ [cm]", Mbin), 540, EtaEdge_min, EtaEdge_max, 450, VtxZEdge_min, VtxZEdge_max)
0380             )
0381         );
0382 
0383         h2D_NClusEtaVtxZ_map.insert( // note : all, outer
0384             std::make_pair(
0385                 Form("h2D_outer_NClusEta_Mbin%d_FineBin", Mbin),
0386                 new TH2D(Form("h2D_outer_NClusEta_Mbin%d_FineBin", Mbin), Form("h2D_outer_NClusEta_Mbin%d_FineBin; Clus #eta (outer);INTT vtxZ [cm]", Mbin), 540, EtaEdge_min, EtaEdge_max, 450, VtxZEdge_min, VtxZEdge_max)
0387             )
0388         );
0389 
0390         h2D_NClusEtaVtxZ_map.insert( // note : typeA, inner
0391             std::make_pair(
0392                 Form("h2D_typeA_inner_NClusEta_Mbin%d_FineBin", Mbin),
0393                 new TH2D(Form("h2D_typeA_inner_NClusEta_Mbin%d_FineBin", Mbin), Form("h2D_typeA_inner_NClusEta_Mbin%d_FineBin; Clus #eta (inner, type A);INTT vtxZ [cm]", Mbin), 540, EtaEdge_min, EtaEdge_max, 450, VtxZEdge_min, VtxZEdge_max)
0394             )
0395         );
0396 
0397         h2D_NClusEtaVtxZ_map.insert( // note : typeA, outer
0398             std::make_pair(
0399                 Form("h2D_typeA_outer_NClusEta_Mbin%d_FineBin", Mbin),
0400                 new TH2D(Form("h2D_typeA_outer_NClusEta_Mbin%d_FineBin", Mbin), Form("h2D_typeA_outer_NClusEta_Mbin%d_FineBin; Clus #eta (outer, type A);INTT vtxZ [cm]", Mbin), 540, EtaEdge_min, EtaEdge_max, 450, VtxZEdge_min, VtxZEdge_max)
0401             )
0402         );
0403     } // note : end of the Mbin loop for h2D    
0404     
0405     h2D_RecoEvtCount_vtxZCentrality = 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%
0406 
0407     h2D_InttVtxZ_Centrality = new TH2D(Form("h2D_InttVtxZ_Centrality"), Form("h2D_InttVtxZ_Centrality;INTT vtxZ [cm];Centrality"), nVtxZBin, VtxZEdge_min, VtxZEdge_max, nCentralityFineBin, CentralityFineEdge_min, CentralityFineEdge_max);
0408 
0409     // Division: ---For Truth-----------------------------------------------------------------------------------------------------------------------------------------
0410     if (runnumber == -1){
0411         h1D_TrueEta_map.clear();
0412         h2D_TrueEtaVtxZ_map.clear();
0413 
0414         for (int Mbin = 0; Mbin < nCentrality_bin; Mbin++)
0415         {
0416             for (int VtxZ_bin = 0; VtxZ_bin < nVtxZBin; VtxZ_bin++)
0417             {
0418                 h1D_TrueEta_map.insert( std::make_pair(
0419                         Form("h1D_TrueEta_Mbin%d_VtxZ%d", Mbin, VtxZ_bin),
0420                         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)
0421                     )
0422                 );
0423             }
0424 
0425             h2D_TrueEtaVtxZ_map.insert( std::make_pair(
0426                     Form("h2D_TrueEtaVtxZ_Mbin%d", Mbin),
0427                     new TH2D(Form("h2D_TrueEtaVtxZ_Mbin%d", Mbin), Form("h2D_TrueEtaVtxZ_Mbin%d;PHG4Particle #eta;TruthPV_trig_z [cm]", Mbin), nEtaBin, EtaEdge_min, EtaEdge_max, nVtxZBin, VtxZEdge_min, VtxZEdge_max)
0428                 )
0429             );
0430 
0431             h2D_TrueEtaVtxZ_map.insert( std::make_pair(
0432                     Form("h2D_TrueEtaVtxZ_Mbin%d_FineBin", Mbin),
0433                     new TH2D(Form("h2D_TrueEtaVtxZ_Mbin%d_FineBin", Mbin), Form("h2D_TrueEtaVtxZ_Mbin%d_FineBin;PHG4Particle #eta;TruthPV_trig_z [cm]", Mbin), 540, EtaEdge_min, EtaEdge_max, 450, VtxZEdge_min, VtxZEdge_max)
0434                 )
0435             );
0436         }
0437 
0438         h2D_TrueEvtCount_vtxZCentrality = new TH2D(Form("h2D_TrueEvtCount_vtxZCentrality"), Form("h2D_TrueEvtCount_vtxZCentrality;TruthPV_trig_z [cm];Centrality [%]"), nVtxZBin, VtxZEdge_min, VtxZEdge_max, nCentrality_bin, &centrality_edges[0]); // note : 0 - 5%
0439     }
0440         
0441 }
0442 
0443 void ClusHistogram::EvtCleanUp()
0444 {
0445     evt_sPH_inner_nocolumn_vec.clear();
0446     evt_sPH_outer_nocolumn_vec.clear();
0447 }
0448 
0449 void ClusHistogram::PrepareClusterVec()
0450 {
0451     if (geo_offset_map.size() == 0){
0452         std::cout<<"the set Geo Offset is set but no input"<<std::endl;
0453         exit(1);
0454     }
0455 
0456     for (int clu_i = 0; clu_i < ClusX -> size(); clu_i++)
0457     {
0458         ClusHistogram::clu_info this_clu;
0459 
0460         this_clu.adc = ClusAdc -> at(clu_i);
0461         this_clu.phi_size = ClusPhiSize -> at(clu_i);
0462         this_clu.sensorZID = ClusLadderZId -> at(clu_i);
0463         this_clu.ladderPhiID = ClusLadderPhiId -> at(clu_i);
0464         this_clu.layerID = ClusLayer -> at(clu_i);
0465 
0466         this_clu.index = clu_i;
0467 
0468         this_clu.x = ClusX -> at(clu_i) + geo_offset_map[Form("%i_%i",this_clu.layerID,this_clu.ladderPhiID)][0];
0469         this_clu.y = ClusY -> at(clu_i) + geo_offset_map[Form("%i_%i",this_clu.layerID,this_clu.ladderPhiID)][1];
0470         this_clu.z = ClusZ -> at(clu_i) + geo_offset_map[Form("%i_%i",this_clu.layerID,this_clu.ladderPhiID)][2];
0471         this_clu.eta_INTTz = get_clu_eta({vertexXYIncm.first, vertexXYIncm.second, INTTvtxZ}, {this_clu.x, this_clu.y, this_clu.z});
0472 
0473         // this_clu.eta_INTTz = ClusEta_INTTz -> at(clu_i);
0474         // this_clu.eta_MBDz  = ClusEta_MBDz  -> at(clu_i);
0475         // this_clu.phi_avgXY = ClusPhi_AvgPV -> at(clu_i);
0476 
0477         // if (runnumber == -1){
0478         //     this_clu.eta_TrueXYZ = ClusEta_TrueXYZ -> at(clu_i);
0479         //     this_clu.phi_TrueXY  = ClusPhi_TrueXY  -> at(clu_i);
0480         // }
0481         // todo: test
0482         // if (this_clu.sensorZID == 0 && (this_clu.layerID == 4) && this_clu.ladderPhiID == 0) {continue;}
0483         // if (this_clu.sensorZID == 0 && (this_clu.layerID == 6) && this_clu.ladderPhiID == 0) {continue;}
0484         // if (this_clu.sensorZID == 0) {continue;}
0485 
0486         if (ColMulMask){
0487             
0488             // note : start from 1 or -1 for outliers 
0489             int GoodColMap_ZId = h1D_GoodColMap_ZId -> Fill(ClusZ -> at(clu_i)); // note : to be free from the geometry offset
0490             if (GoodColMap_ZId == -1) {continue;}
0491 
0492             int GoodColMap_XId = (this_clu.layerID - 3) * 20 + this_clu.ladderPhiID + 1; // note : Layer 4, phiID 9 -> 30
0493 
0494             if (
0495                 h2D_GoodColMap != nullptr && 
0496                 h2D_GoodColMap -> GetBinContent(GoodColMap_XId, GoodColMap_ZId) == 0
0497             )
0498             {
0499                 continue;
0500             }
0501         }
0502 
0503         if (isClusQA.first && this_clu.adc <= isClusQA.second.first) {continue;} // note : adc
0504         if (isClusQA.first && this_clu.phi_size > isClusQA.second.second) {continue;} // note : phi size
0505 
0506         std::vector<ClusHistogram::clu_info>* p_evt_sPH_nocolumn_vec =
0507         (this_clu.layerID == 3 || this_clu.layerID == 4) ? (&evt_sPH_inner_nocolumn_vec) : (&evt_sPH_outer_nocolumn_vec);
0508         p_evt_sPH_nocolumn_vec -> push_back(this_clu);
0509     }
0510 
0511     if (SetRandomHits.first){
0512         if (inner_UniqueClusXYZ_map.size() == 0 || outer_UniqueClusXYZ_map.size() == 0){
0513             std::cout<<"the Unique Clus XYZ is not set for generating the random hits"<<std::endl;
0514             exit(1);
0515         }
0516 
0517         for (int clu_i = 0; clu_i < SetRandomHits.second; clu_i++)
0518         {
0519             ClusHistogram::clu_info this_clu;
0520 
0521             int inner_or_outer = (rand() % 2 == 0) ? 0 : 1; // note : 0 -> inner, 1 -> outer
0522 
0523             std::map<std::string, std::tuple<double, double, double, int, int, int, int, int>>* p_UniqueClusXYZ_map = (inner_or_outer == 0) ? (&inner_UniqueClusXYZ_map) : (&outer_UniqueClusXYZ_map);
0524             std::vector<std::string>* p_UniqueClusXYZ_vec                                                           = (inner_or_outer == 0) ? (&inner_UniqueClusXYZ_vec) : (&outer_UniqueClusXYZ_vec);
0525             int Rand_index = int(rand3 -> Uniform(0, p_UniqueClusXYZ_vec -> size()));
0526 
0527             double selected_x        = std::get<0>(p_UniqueClusXYZ_map -> at(p_UniqueClusXYZ_vec -> at(Rand_index)));
0528             double selected_y        = std::get<1>(p_UniqueClusXYZ_map -> at(p_UniqueClusXYZ_vec -> at(Rand_index)));
0529             double selected_z        = std::get<2>(p_UniqueClusXYZ_map -> at(p_UniqueClusXYZ_vec -> at(Rand_index)));
0530             int selected_sensorZID   = std::get<3>(p_UniqueClusXYZ_map -> at(p_UniqueClusXYZ_vec -> at(Rand_index)));
0531             int selected_layerID     = std::get<4>(p_UniqueClusXYZ_map -> at(p_UniqueClusXYZ_vec -> at(Rand_index)));
0532             int selected_adc         = std::get<5>(p_UniqueClusXYZ_map -> at(p_UniqueClusXYZ_vec -> at(Rand_index)));
0533             int selected_phi_size    = std::get<6>(p_UniqueClusXYZ_map -> at(p_UniqueClusXYZ_vec -> at(Rand_index)));
0534             int selected_ladderPhiID = std::get<7>(p_UniqueClusXYZ_map -> at(p_UniqueClusXYZ_vec -> at(Rand_index)));
0535 
0536             this_clu.x = selected_x;
0537             this_clu.y = selected_y;
0538             this_clu.z = selected_z;
0539             this_clu.eta_INTTz = get_clu_eta({vertexXYIncm.first, vertexXYIncm.second, INTTvtxZ}, {selected_x, selected_y, selected_z});
0540             this_clu.sensorZID = selected_sensorZID;
0541             this_clu.layerID = selected_layerID;
0542             this_clu.adc = selected_adc;
0543             this_clu.phi_size = selected_phi_size;
0544             this_clu.ladderPhiID = selected_ladderPhiID;
0545 
0546             this_clu.index = ClusX -> size() + clu_i;
0547 
0548             std::vector<ClusHistogram::clu_info>* p_evt_sPH_nocolumn_vec =
0549             (inner_or_outer == 0) ? (&evt_sPH_inner_nocolumn_vec) : (&evt_sPH_outer_nocolumn_vec);
0550             p_evt_sPH_nocolumn_vec -> push_back(this_clu);
0551         }
0552     }
0553 }
0554 
0555 void ClusHistogram::PrepareUniqueClusXYZ()
0556 {
0557     std::cout<<"In preparaing the unique cluster XYZ map"<<std::endl;
0558 
0559     inner_UniqueClusXYZ_map.clear();
0560     outer_UniqueClusXYZ_map.clear();
0561     inner_UniqueClusXYZ_vec.clear();
0562     outer_UniqueClusXYZ_vec.clear();
0563 
0564     h2D_RandClusXY_ref = new TH2D("h2D_RandClusXY_ref","h2D_RandClusXY_ref;ClusX [cm];ClusY [cm]",200,-15,15,200,-15,15);
0565     h1D_RandClusZ_ref = new TH1D("h1D_RandClusZ_ref","h1D_RandClusZ_ref;ClusZ [cm];Entries",100,-30,30);
0566 
0567     if (geo_offset_map.size() == 0){
0568         std::cout<<"the set Geo Offset is set but no input"<<std::endl;
0569         exit(1);
0570     }
0571 
0572     for (int i = 0; i < 500; i++) // todo : but it's limited by the MC sample
0573     {
0574         tree_in -> GetEntry(i);
0575 
0576         // std::cout<<111<<std::endl;
0577 
0578         for (int clu_i = 0; clu_i < ClusZ -> size(); clu_i++){
0579             
0580             // std::cout<<222<<std::endl;
0581             if (ClusPhiSize -> at(clu_i) > 1) {continue;} // todo: maybe we don't need it
0582             // std::cout<<333<<std::endl;
0583 
0584             double clu_x = ClusX -> at(clu_i) + geo_offset_map[Form("%i_%i",ClusLayer->at(clu_i),ClusLadderPhiId->at(clu_i))][0];
0585             double clu_y = ClusY -> at(clu_i) + geo_offset_map[Form("%i_%i",ClusLayer->at(clu_i),ClusLadderPhiId->at(clu_i))][1];
0586             double clu_z = ClusZ -> at(clu_i) + geo_offset_map[Form("%i_%i",ClusLayer->at(clu_i),ClusLadderPhiId->at(clu_i))][2];
0587             
0588             // std::cout<<4444<<std::endl;
0589 
0590             // std::string clu_key = Form("%.2f_%.2f_%.2f", clu_x, clu_y, clu_z);
0591             std::string clu_key = Form("%.2f_%.2f_%d", clu_x, clu_y, h1D_RandClusZ_ref->FindBin(clu_z));
0592 
0593             // std::cout<<555<<std::endl;
0594 
0595             std::map<std::string, std::tuple<double, double, double, int, int, int, int, int>>* p_UniqueClusXYZ_map = (ClusLayer -> at(clu_i) == 3 || ClusLayer -> at(clu_i) == 4) ? (&inner_UniqueClusXYZ_map) : (&outer_UniqueClusXYZ_map);
0596             std::vector<std::string>* p_UniqueClusXYZ_vec                                                           = (ClusLayer -> at(clu_i) == 3 || ClusLayer -> at(clu_i) == 4) ? (&inner_UniqueClusXYZ_vec) : (&outer_UniqueClusXYZ_vec);
0597 
0598             // std::cout<<666<<std::endl;
0599 
0600             if (ColMulMask){
0601                 
0602                 // note : start from 1 or -1 for outliers 
0603                 int GoodColMap_ZId = h1D_GoodColMap_ZId -> Fill(ClusZ -> at(clu_i)); // note : to be free from the geometry offset
0604                 if (GoodColMap_ZId == -1) {continue;}
0605 
0606                 int GoodColMap_XId = (ClusLayer->at(clu_i) - 3) * 20 + ClusLadderPhiId->at(clu_i) + 1; // note : Layer 4, phiID 9 -> 30
0607 
0608                 if (
0609                     h2D_GoodColMap != nullptr && 
0610                     h2D_GoodColMap -> GetBinContent(GoodColMap_XId, GoodColMap_ZId) == 0
0611                 )
0612                 {
0613                     continue;
0614                 }
0615             }
0616 
0617             if (p_UniqueClusXYZ_map -> find(clu_key) == p_UniqueClusXYZ_map -> end()){
0618                 p_UniqueClusXYZ_map -> insert( 
0619                     std::make_pair(
0620                         clu_key,
0621                         std::make_tuple(
0622                             clu_x, clu_y, clu_z, 
0623                             ClusLadderZId->at(clu_i), 
0624                             ClusLayer -> at(clu_i),
0625                             ClusAdc -> at(clu_i),
0626                             ClusPhiSize -> at(clu_i),
0627                             ClusLadderPhiId -> at(clu_i)
0628                         )
0629                     ) 
0630                 );
0631                 
0632                 // std::cout<<7777<<std::endl;
0633 
0634                 p_UniqueClusXYZ_vec -> push_back(clu_key);
0635 
0636                 // std::cout<<8888<<std::endl;
0637 
0638                 h2D_RandClusXY_ref -> Fill(clu_x, clu_y);
0639                 h1D_RandClusZ_ref -> Fill(clu_z);
0640 
0641                 // std::cout<<9999<<std::endl;
0642             }            
0643         }
0644     }
0645 }
0646 
0647 double ClusHistogram::CheckGeoOffsetMap()
0648 {
0649     double sum = 0;
0650     for (auto &pair : geo_offset_map)
0651     {
0652         sum += fabs(pair.second[0]) + fabs(pair.second[1]) + fabs(pair.second[2]);
0653     }
0654     return sum;
0655 }
0656 
0657 void ClusHistogram::MainProcess()
0658 {
0659     if (SetRandomHits.first){PrepareUniqueClusXYZ();}
0660 
0661     if (ColMulMask && h2D_GoodColMap == nullptr){
0662         std::cout<<"The GoodColMap is not set correctly"<<std::endl;
0663         exit(1);
0664     }
0665 
0666     if (
0667         ColMulMask &&
0668         (
0669             h1D_GoodColMap_ZId -> GetNbinsX() != h2D_GoodColMap -> GetNbinsY() ||
0670             fabs(h1D_GoodColMap_ZId -> GetXaxis() -> GetXmin() - h2D_GoodColMap -> GetYaxis() -> GetXmin()) > 0.0000001 ||
0671             fabs(h1D_GoodColMap_ZId -> GetXaxis() -> GetXmax() - h2D_GoodColMap -> GetYaxis() -> GetXmax()) > 0.0000001
0672         )
0673 
0674     ){
0675         std::cout<<"The setting of h1D_GoodColMap_ZId is different from h2D_GoodColMap"<<std::endl;
0676         std::cout<<"h1D_GoodColMap_ZId : "<<h1D_GoodColMap_ZId -> GetNbinsX()<<" "<<h1D_GoodColMap_ZId -> GetXaxis() -> GetXmin()<<" "<<h1D_GoodColMap_ZId -> GetXaxis() -> GetXmax()<<std::endl;
0677         std::cout<<"h2D_GoodColMap : "<<h2D_GoodColMap -> GetNbinsY()<<" "<<h2D_GoodColMap -> GetYaxis() -> GetXmin()<<" "<<h2D_GoodColMap -> GetYaxis() -> GetXmax()<<std::endl;
0678         exit(1);
0679     }
0680 
0681     if (HaveGeoOffsetTag && CheckGeoOffsetMap() <= 0){
0682         std::cout<<"The geo offset map is not set correctly"<<std::endl;
0683         exit(1);
0684     }
0685 
0686     if (HaveGeoOffsetTag == false && CheckGeoOffsetMap() > 0.0001) {
0687         std::cout<<"The HaveGeoOffsetTag is false, but the GeoOffsetMap has some non-zero numbers, please check the GeoOffsetMap"<<std::endl;
0688         exit(1);
0689     }
0690 
0691     for (int i = 0; i < run_nEvents; i++)
0692     {
0693         tree_in -> GetEntry(i);
0694 
0695         EvtCleanUp();
0696 
0697         if (RandInttZ){
0698             INTTvtxZ = rand3 -> Uniform(VtxZEdge_min,VtxZEdge_max);
0699             INTTvtxZError = 0;
0700         }
0701 
0702         if (i % 10 == 0) {std::cout << "Processing event " << i<<", NClus : "<< ClusX -> size() << std::endl;}
0703 
0704         // =======================================================================================================================================================
0705         // note : hard cut
0706 
0707         // note : for data
0708         if (runnumber != -1 && BcoFullDiffCut && InttBcoFullDiff_next <= cut_InttBcoFullDIff_next) {continue;}
0709         if (runnumber != -1 && MBDNSg2 != 1) {continue;} // todo: assume MC no trigger
0710 
0711         // note : for MC
0712         // if (runnumber == -1 && NTruthVtx != 1) {continue;}
0713 
0714         // note : both data and MC
0715         if (is_min_bias != 1) {continue;}
0716         if (MBD_z_vtx != MBD_z_vtx) {continue;}
0717         if (MBD_centrality != MBD_centrality) {continue;}
0718         if (MBD_centrality < 0 || MBD_centrality > 100) {continue;}
0719         if (INTTvtxZ != INTTvtxZ) {continue;}
0720         if (MBD_z_vtx < cut_GlobalMBDvtxZ.first || MBD_z_vtx > cut_GlobalMBDvtxZ.second) {continue;} // todo: the hard cut 60 cm 
0721         // =======================================================================================================================================================
0722 
0723         int Mbin = h1D_centrality_bin -> Fill(MBD_centrality);
0724         Mbin = (Mbin == -1) ? -1 : Mbin - 1;
0725         if (Mbin == -1) {
0726             std::cout << "Mbin == -1, MBD_centrality = " << MBD_centrality << std::endl;
0727             continue;
0728         }
0729         // =======================================================================================================================================================
0730 
0731 
0732         // note : for truth
0733         if (runnumber == -1){
0734             int TruthVtxZ_bin = h1D_vtxz_template->Fill(TruthPV_trig_z);
0735             TruthVtxZ_bin = (TruthVtxZ_bin == -1) ? -1 : TruthVtxZ_bin - 1;
0736 
0737             if (TruthVtxZ_bin != -1){
0738                 for (int true_i = 0; true_i < NPrimaryG4P; true_i++){
0739                     if (PrimaryG4P_isChargeHadron->at(true_i) != 1) { continue; }
0740                     
0741                     h1D_TrueEta_map[Form("h1D_TrueEta_Mbin%d_VtxZ%d", Mbin, TruthVtxZ_bin)] -> Fill(PrimaryG4P_Eta->at(true_i));
0742                     h2D_TrueEtaVtxZ_map[Form("h2D_TrueEtaVtxZ_Mbin%d", Mbin)] -> Fill(PrimaryG4P_Eta->at(true_i), TruthPV_trig_z);
0743                     h2D_TrueEtaVtxZ_map[Form("h2D_TrueEtaVtxZ_Mbin%d_FineBin", Mbin)] -> Fill(PrimaryG4P_Eta->at(true_i), TruthPV_trig_z);
0744                 }
0745 
0746                 h2D_TrueEvtCount_vtxZCentrality->Fill(TruthPV_trig_z, MBD_centrality);
0747             }
0748         }
0749 
0750 
0751         // =======================================================================================================================================================
0752         // note : optional cut
0753         if (INTT_vtxZ_QA && (MBD_z_vtx - INTTvtxZ < cut_vtxZDiff.first || MBD_z_vtx - INTTvtxZ > cut_vtxZDiff.second) ) {continue;}
0754         if (INTT_vtxZ_QA && (TrapezoidalFitWidth < cut_TrapezoidalFitWidth.first || TrapezoidalFitWidth > cut_TrapezoidalFitWidth.second)){continue;}
0755         if (INTT_vtxZ_QA && (TrapezoidalFWHM < cut_TrapezoidalFWHM.first || TrapezoidalFWHM > cut_TrapezoidalFWHM.second)){continue;}
0756         if (INTT_vtxZ_QA && (INTTvtxZError < cut_INTTvtxZError.first || INTTvtxZError > cut_INTTvtxZError.second)){continue;}
0757         // =======================================================================================================================================================
0758         
0759         if (vtxZReweight.first && runnumber != -1){
0760             std::cout<<"Should not have the vtxZ weighting from the data"<<std::endl;
0761             exit(1);
0762         }
0763 
0764         double INTTvtxZWeighting;
0765         if (vtxZReweight.first && h1D_INTT_vtxZ_reweighting != nullptr){
0766             INTTvtxZWeighting = h1D_INTT_vtxZ_reweighting -> GetBinContent(h1D_INTT_vtxZ_reweighting -> FindBin(INTTvtxZ));
0767         }
0768         else if (vtxZReweight.first && h1D_INTT_vtxZ_reweighting == nullptr){
0769             std::cout << "ApplyVtxZReWeighting is true, but h1D_INTT_vtxZ_reweighting is nullptr" << std::endl;
0770             exit(1);
0771         }
0772         else {
0773             INTTvtxZWeighting = 1.0;
0774         }
0775 
0776         // =======================================================================================================================================================
0777         
0778         int InttVtxZ_bin = h1D_vtxz_template->Fill(INTTvtxZ);
0779         InttVtxZ_bin = (InttVtxZ_bin == -1) ? -1 : InttVtxZ_bin - 1;
0780 
0781         // int MBDVtxZ_bin = h1D_vtxz_template->Fill(MBD_z_vtx);
0782         // MBDVtxZ_bin = (MBDVtxZ_bin == -1) ? -1 : MBDVtxZ_bin - 1;
0783 
0784         // note : everything below is within INTT vtxZ -45 cm  ~ 45 cm
0785         if (InttVtxZ_bin == -1) {continue;}
0786 
0787 
0788         // note : for reconstructed 
0789         // note : INTT vtxZ range -45 cm ~ 45 cm
0790         h1D_centrality_bin_weighted -> Fill(MBD_centrality,INTTvtxZWeighting);
0791         h1D_INTTvtxZ -> Fill(INTTvtxZ, INTTvtxZWeighting);
0792         h1D_centrality_map["h1D_centrality"] -> Fill(MBD_centrality,INTTvtxZWeighting);
0793 
0794         h2D_RecoEvtCount_vtxZCentrality -> Fill(INTTvtxZ, MBD_centrality, INTTvtxZWeighting);
0795         h2D_InttVtxZ_Centrality -> Fill(INTTvtxZ, MBD_centrality, INTTvtxZWeighting); // note : fine bins in centrality
0796         h1D_centrality_map[Form("h1D_centrality_InttVtxZ%d", InttVtxZ_bin)] -> Fill(MBD_centrality,INTTvtxZWeighting);
0797 
0798         PrepareClusterVec();
0799 
0800         // note : inner 
0801         for (int clu_i = 0; clu_i < evt_sPH_inner_nocolumn_vec.size(); clu_i++){
0802             ClusHistogram::clu_info this_clu = evt_sPH_inner_nocolumn_vec[clu_i];
0803 
0804             h1D_NClusEta_map[Form("h1D_inner_NClusEta_Mbin%d_VtxZ%d", Mbin, InttVtxZ_bin)] -> Fill(this_clu.eta_INTTz, INTTvtxZWeighting);
0805             h2D_NClusEtaVtxZ_map[Form("h2D_inner_NClusEta_Mbin%d", Mbin)] -> Fill(this_clu.eta_INTTz, INTTvtxZ, INTTvtxZWeighting);
0806             h2D_NClusEtaVtxZ_map[Form("h2D_inner_NClusEta_Mbin%d_FineBin", Mbin)] -> Fill(this_clu.eta_INTTz, INTTvtxZ, INTTvtxZWeighting);
0807             
0808             // note : type A
0809             if (this_clu.sensorZID == typeA_sensorZID[0] || this_clu.sensorZID == typeA_sensorZID[1]){
0810                 h1D_NClusEta_map[Form("h1D_typeA_inner_NClusEta_Mbin%d_VtxZ%d", Mbin, InttVtxZ_bin)] -> Fill(this_clu.eta_INTTz, INTTvtxZWeighting);
0811                 h2D_NClusEtaVtxZ_map[Form("h2D_typeA_inner_NClusEta_Mbin%d", Mbin)] -> Fill(this_clu.eta_INTTz, INTTvtxZ, INTTvtxZWeighting);
0812                 h2D_NClusEtaVtxZ_map[Form("h2D_typeA_inner_NClusEta_Mbin%d_FineBin", Mbin)] -> Fill(this_clu.eta_INTTz, INTTvtxZ, INTTvtxZWeighting);
0813             }            
0814         }
0815 
0816         // note : outer 
0817         for (int clu_i = 0; clu_i < evt_sPH_outer_nocolumn_vec.size(); clu_i++){
0818             ClusHistogram::clu_info this_clu = evt_sPH_outer_nocolumn_vec[clu_i];
0819 
0820             h1D_NClusEta_map[Form("h1D_outer_NClusEta_Mbin%d_VtxZ%d", Mbin, InttVtxZ_bin)] -> Fill(this_clu.eta_INTTz, INTTvtxZWeighting);
0821             h2D_NClusEtaVtxZ_map[Form("h2D_outer_NClusEta_Mbin%d", Mbin)] -> Fill(this_clu.eta_INTTz, INTTvtxZ, INTTvtxZWeighting);
0822             h2D_NClusEtaVtxZ_map[Form("h2D_outer_NClusEta_Mbin%d_FineBin", Mbin)] -> Fill(this_clu.eta_INTTz, INTTvtxZ, INTTvtxZWeighting);
0823             
0824             // note : type A
0825             if (this_clu.sensorZID == typeA_sensorZID[0] || this_clu.sensorZID == typeA_sensorZID[1]){
0826                 h1D_NClusEta_map[Form("h1D_typeA_outer_NClusEta_Mbin%d_VtxZ%d", Mbin, InttVtxZ_bin)] -> Fill(this_clu.eta_INTTz, INTTvtxZWeighting);
0827                 h2D_NClusEtaVtxZ_map[Form("h2D_typeA_outer_NClusEta_Mbin%d", Mbin)] -> Fill(this_clu.eta_INTTz, INTTvtxZ, INTTvtxZWeighting);
0828                 h2D_NClusEtaVtxZ_map[Form("h2D_typeA_outer_NClusEta_Mbin%d_FineBin", Mbin)] -> Fill(this_clu.eta_INTTz, INTTvtxZ, INTTvtxZWeighting);
0829             }
0830             
0831         }
0832 
0833 
0834     } // note : end pf event looping
0835 }
0836 
0837 double ClusHistogram::get_clu_eta(std::vector<double> vertex, std::vector<double> clu_pos)
0838 {
0839     double correct_x = clu_pos[0] - vertex[0];
0840     double correct_y = clu_pos[1] - vertex[1];
0841     double correct_z = clu_pos[2] - vertex[2];
0842     double clu_r = sqrt(pow(correct_x,2) + pow(correct_y,2));
0843     // cout<<"correct info : "<<correct_x<<" "<<correct_y<<" "<<correct_z<<" "<<clu_r<<endl;
0844 
0845     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)));
0846 }
0847 
0848 void ClusHistogram::EndRun()
0849 {
0850     file_out -> cd();
0851     tree_out_par -> Fill();
0852     tree_out_par -> Write();
0853 
0854     if (h1D_INTT_vtxZ_reweighting != nullptr){
0855         h1D_INTT_vtxZ_reweighting -> Write("h1D_Used_INTTvtxZ_reweighting");
0856     }
0857 
0858     h1D_vtxz_template -> Write();
0859     h1D_INTTvtxZ -> Write();
0860     h1D_centrality_bin -> Write();
0861     h1D_centrality_bin_weighted -> Write();
0862     h2D_RecoEvtCount_vtxZCentrality -> Write();
0863     h2D_InttVtxZ_Centrality -> Write();
0864     if (runnumber == -1) {h2D_TrueEvtCount_vtxZCentrality->Write();}
0865 
0866 
0867     for (auto &pair : h2D_NClusEtaVtxZ_map){
0868         pair.second -> Write();
0869     }
0870 
0871 
0872     if (runnumber == -1){
0873         for (auto &pair : h2D_TrueEtaVtxZ_map){
0874             pair.second -> Write();
0875         }
0876     }
0877 
0878 
0879     for (auto &pair : h1D_NClusEta_map){
0880         pair.second -> Write();
0881     }
0882 
0883 
0884     for (auto &pair : h1D_centrality_map){
0885         pair.second -> Write();
0886     }
0887 
0888     
0889     if (runnumber == -1){
0890         for (auto &pair : h1D_TrueEta_map){
0891             pair.second -> Write();
0892         }
0893     }
0894 
0895     if (h1D_RandClusZ_ref != nullptr && h2D_RandClusXY_ref != nullptr){
0896         h2D_RandClusXY_ref -> Write();
0897         h1D_RandClusZ_ref -> Write();
0898     }
0899 
0900     
0901     file_out -> Close();
0902 }
[ 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!