Back to home page

sPhenix code displayed by LXR

 
 

    


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

0001 #include "ColumnCheck.h"
0002 
0003 ColumnCheck::ColumnCheck(
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     double SetMbinFloat_in,
0014 
0015     std::pair<double, double> VtxZRange_in,
0016     bool IsZClustering_in,
0017     bool BcoFullDiffCut_in,
0018     std::pair<bool, std::pair<double, double>> isClusQA_in, // note : {adc, phi size}
0019 
0020     bool ColMulMask_in
0021 
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     {false, nullptr}, // note : vtxZReweight_in
0034     BcoFullDiffCut_in,
0035     false, // note : INTT_vtxZ_QA_in
0036     isClusQA_in, // note : {adc, phi size}
0037     false, // note : HaveGeoOffsetTag_in
0038     {false, 0}, // note : SetRandomHits_in
0039     false, // note : RandInttZ_in
0040     ColMulMask_in,
0041 
0042     1 // note : c_type_in
0043 ), 
0044 IsZClustering(IsZClustering_in),
0045 VtxZRange(VtxZRange_in),
0046 SetMbinFloat(SetMbinFloat_in) // note : 0-100
0047 
0048 {
0049     PrepareOutPutFileName();
0050     PrepareOutPutRootFile();
0051     PrepareHistograms();
0052 }
0053 
0054 void ColumnCheck::PrepareOutPutFileName()
0055 {
0056     std::string job_index = std::to_string( process_id );
0057     int job_index_len = 5;
0058     job_index.insert(0, job_index_len - job_index.size(), '0');
0059 
0060     std::string runnumber_str = std::to_string( runnumber );
0061     if (runnumber != -1){
0062         int runnumber_str_len = 8;
0063         runnumber_str.insert(0, runnumber_str_len - runnumber_str.size(), '0');
0064     }
0065 
0066     if (output_file_name_suffix.size() > 0 && output_file_name_suffix[0] != '_') {
0067         output_file_name_suffix = "_" + output_file_name_suffix;
0068     }
0069 
0070     output_filename = "ColumnCheck";
0071     output_filename = (runnumber != -1) ? "Data_" + output_filename : "MC_" + output_filename;
0072     output_filename += (BcoFullDiffCut && runnumber != -1) ? "_BcoFullDiffCut" : "";
0073     output_filename += (IsZClustering) ? "_ZClustering" : "";
0074     output_filename += Form("_Mbin%.0f",SetMbinFloat);
0075     output_filename += (ColMulMask) ? "_ColMulMask" : "";
0076     output_filename += Form("_VtxZ%.0fto%.0fcm",VtxZRange.first,VtxZRange.second);
0077     output_filename += (isClusQA.first) ? Form("_ClusQAAdc%.0fPhiSize%.0f",isClusQA.second.first,isClusQA.second.second) : "";
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 ColumnCheck::PrepareOutPutRootFile()
0083 {
0084     file_out = new TFile(Form("%s/%s", output_directory.c_str(), output_filename.c_str()), "RECREATE");
0085 }
0086 
0087 void ColumnCheck::PrepareHistograms()
0088 {
0089     h1D_map.clear();
0090     h2D_map.clear();
0091 
0092     h1D_GoodColMap_ZId = new TH1D("h1D_GoodColMap_ZId","h1D_GoodColMap_ZId;ClusZ [cm];Entries",nZbin, Zmin, Zmax);
0093 
0094     h1D_map["h1D_ClusZ"] = new TH1D("h1D_ClusZ","h1D_ClusZ;ClusZ [cm];Entries",nZbin, Zmin, Zmax);
0095 
0096     for (int i = 0; i < nZbin; i++)
0097     {
0098         h2D_map[Form("h2D_ClusCountLayerPhiId_ZId%d",i)] = new TH2D(Form("h2D_ClusCountLayerPhiId_ZId%d",i),Form("h2D_ClusCountLayerPhiId_ZId%d;LayerID;LadderPhiId",i),4,3,7,16,0,16);
0099     }
0100 
0101     for (int layer_i = 3; layer_i < 7; layer_i++)
0102     {   
0103         int phi_max = (layer_i == 3 || layer_i == 4) ? 12 : 16;
0104 
0105         for (int phi_i = 0; phi_i < phi_max; phi_i++)
0106         {
0107             h2D_map[Form("h2D_ClusRadiusZID_Layer%d_PhiId%d",layer_i,phi_i)] = new TH2D(Form("h2D_ClusRadiusZID_Layer%d_PhiId%d",layer_i,phi_i),Form("h2D_ClusRadiusZID_Layer%d_PhiId%d;Clus radius [cm]; ClusZ [cm]",layer_i,phi_i), 275, 6, 11.5, nZbin, Zmin, Zmax);
0108         }
0109     }
0110 
0111     h2D_map[Form("h2D_ClusRadiusZID_All")] = new TH2D(Form("h2D_ClusRadiusZID_All"),Form("h2D_ClusRadiusZID_All;Clus radius [cm]; ClusZ [cm]"), 275, 6, 11.5, nZbin, Zmin, Zmax);
0112 
0113     for (int i = 0; i < 4; i++){ // note : layer
0114         h2D_map[Form("h2D_ClusCountPhiIdZId_Layer%d",i+3)] = new TH2D(Form("h2D_ClusCountPhiIdZId_Layer%d",i+3),Form("h2D_ClusCountPhiIdZId_Layer%d;LadderPhiId;ZId",i+3),16,0,16, nZbin, Zmin, Zmax);
0115     }
0116 }
0117 
0118 void ColumnCheck::EvtCleanUp()
0119 {
0120     evt_sPH_inner_nocolumn_vec.clear();
0121     evt_sPH_outer_nocolumn_vec.clear();   
0122 }
0123 
0124 void ColumnCheck::MainProcess()
0125 {
0126     for (int i = 0; i < run_nEvents; i++)
0127     {
0128         tree_in -> GetEntry(i);
0129 
0130         EvtCleanUp();
0131 
0132         if (i % 10 == 0) {std::cout << "Processing event " << i<<", NClus : "<< ClusX -> size() << std::endl;}
0133 
0134         // =======================================================================================================================================================
0135         // note : hard cut
0136 
0137         // note : for data
0138         if (runnumber != -1 && BcoFullDiffCut && InttBcoFullDiff_next <= cut_InttBcoFullDIff_next) {continue;}
0139         if (runnumber != -1 && MBDNSg2 != 1) {continue;} // todo: assume MC no trigger
0140 
0141         // note : for MC
0142         // if (runnumber == -1 && NTruthVtx != 1) {continue;}
0143 
0144         // note : both data and MC
0145         if (is_min_bias != 1) {continue;}
0146         if (MBD_z_vtx != MBD_z_vtx) {continue;}
0147         if (MBD_centrality != MBD_centrality) {continue;}
0148         if (MBD_centrality < 0 || MBD_centrality > 100) {continue;}
0149         // if (INTTvtxZ != INTTvtxZ) {continue;}
0150 
0151         if (MBD_centrality > SetMbinFloat) {continue;}
0152         if (MBD_z_vtx < VtxZRange.first || MBD_z_vtx > VtxZRange.second) {continue;}
0153         // =======================================================================================================================================================
0154 
0155         PrepareClusterVec();
0156 
0157         for (int barrel_i = 0; barrel_i < 2; barrel_i++)
0158         {
0159             std::vector<ClusHistogram::clu_info> &this_vec = (barrel_i == 0) ? evt_sPH_inner_nocolumn_vec : evt_sPH_outer_nocolumn_vec;
0160 
0161             for (ClusHistogram::clu_info this_clu : this_vec)
0162             {
0163                 int ZID = h1D_map["h1D_ClusZ"] -> Fill(this_clu.z) - 1;
0164                 if (ZID == -2) {continue;}
0165 
0166                 double ClusRadius = sqrt(pow(this_clu.x - vertexXYIncm.first,2) + pow(this_clu.y - vertexXYIncm.second,2));
0167                 double ColumnLength = (this_clu.sensorZID == typeA_sensorZID[0] || this_clu.sensorZID == typeA_sensorZID[1]) ? typeA_sensor_half_length_incm : typeB_sensor_half_length_incm;
0168 
0169                 h2D_map[Form("h2D_ClusCountLayerPhiId_ZId%d",ZID)] -> Fill(this_clu.layerID, this_clu.ladderPhiID, 1. / (ColumnLength));
0170                 h2D_map[Form("h2D_ClusCountPhiIdZId_Layer%d",this_clu.layerID)] -> Fill(this_clu.ladderPhiID, this_clu.z, 1. / (ColumnLength));
0171                 h2D_map[Form("h2D_ClusRadiusZID_Layer%d_PhiId%d",this_clu.layerID,this_clu.ladderPhiID)] -> Fill(ClusRadius, this_clu.z, 1. / (ColumnLength));
0172                 h2D_map[Form("h2D_ClusRadiusZID_All")] -> Fill(ClusRadius, this_clu.z, 1. / (ColumnLength));
0173             }
0174         }
0175     }
0176 }
0177 
0178 void ColumnCheck::EndRun()
0179 {
0180     file_out -> cd();
0181 
0182     for (auto &pair : h1D_map)
0183     {
0184         pair.second -> Write();
0185     }
0186 
0187     for (auto &pair : h2D_map)
0188     {
0189         pair.second -> Write();
0190     }
0191 
0192     file_out -> Close();
0193 
0194 }