Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "SemiChipClustering.h"
0002 
0003 SemiChipClustering::SemiChipClustering(
0004     int process_id_in, // note : 1 or 2
0005     int runnumber_in, // note : still, (54280 or -1)
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 
0013     bool BcoFullDiffCut_in,
0014     bool INTT_vtxZ_QA_in,
0015     std::pair<bool, std::string> BadChMask_in,
0016     bool ApplyHitQA_in,
0017     bool clone_hit_remove_BCO_tag_in,
0018     std::pair<bool, int> cut_HitBcoDiff_in,
0019     std::vector<int> adc_conversion_vec_in
0020 ):
0021     process_id(process_id_in),
0022     runnumber(runnumber_in),
0023     run_nEvents(run_nEvents_in),
0024     input_directory(input_directory_in),
0025     input_file_name(input_file_name_in),
0026     output_directory(output_directory_in),
0027     output_file_name_suffix(output_file_name_suffix_in),
0028     BcoFullDiffCut(BcoFullDiffCut_in),
0029     INTT_vtxZ_QA(INTT_vtxZ_QA_in),
0030     BadChMask(BadChMask_in),
0031     ApplyHitQA(ApplyHitQA_in),
0032     clone_hit_remove_BCO_tag(clone_hit_remove_BCO_tag_in),
0033     cut_HitBcoDiff(cut_HitBcoDiff_in),
0034     adc_conversion_vec(adc_conversion_vec_in)
0035 {
0036     PrepareInputRootFile();
0037 
0038     run_nEvents = (run_nEvents == -1) ? tree_in->GetEntries() : run_nEvents;
0039     run_nEvents = (run_nEvents > tree_in->GetEntries()) ? tree_in->GetEntries() : run_nEvents;
0040 
0041     PrepareOutPutFileName();
0042     PrepareOutPutRootFile();
0043     
0044     hot_channel_map.clear();
0045     if (BadChMask.first){PrepareHotChannel();}
0046 
0047     PrepareHists();
0048 }
0049 
0050 std::map<std::string, int> SemiChipClustering::GetInputTreeBranchesMap(TTree * m_tree_in)
0051 {
0052     std::map<std::string, int> branch_map;
0053     TObjArray * branch_list = m_tree_in -> GetListOfBranches();
0054     for (int i = 0; i < branch_list -> GetEntries(); i++)
0055     {
0056         TBranch * branch = dynamic_cast<TBranch*>(branch_list->At(i));
0057         branch_map[branch -> GetName()] = 1;
0058     }
0059     return branch_map;
0060 }
0061 
0062 void SemiChipClustering::PrepareInputRootFile(){
0063     file_in = TFile::Open(Form("%s/%s", input_directory.c_str(), input_file_name.c_str()));
0064     if (!file_in || file_in -> IsZombie() || file_in == nullptr) {
0065         std::cout << "Error: cannot open file: " << input_file_name << std::endl;
0066         exit(1);
0067     }
0068 
0069     tree_in = (TTree*)file_in -> Get(tree_name.c_str());
0070     
0071     std::map<std::string, int> branch_map = GetInputTreeBranchesMap(tree_in);
0072 
0073     tree_in -> SetBranchStatus("*", 0);
0074 
0075     tree_in -> SetBranchStatus("MBD_z_vtx", 1);
0076     tree_in -> SetBranchStatus("is_min_bias", 1);
0077     tree_in -> SetBranchStatus("MBD_centrality", 1);
0078     tree_in -> SetBranchStatus("MBD_south_charge_sum", 1);
0079     tree_in -> SetBranchStatus("MBD_north_charge_sum", 1);
0080     tree_in -> SetBranchStatus("MBD_charge_sum", 1);
0081     tree_in -> SetBranchStatus("MBD_charge_asymm", 1);
0082     tree_in -> SetBranchStatus("InttBcoFullDiff_next", 1);
0083 
0084     tree_in -> SetBranchStatus("MBDNSg2",1);
0085     tree_in -> SetBranchStatus("MBDNSg2_vtxZ10cm",1);
0086     tree_in -> SetBranchStatus("MBDNSg2_vtxZ30cm",1);
0087     tree_in -> SetBranchStatus("MBDNSg2_vtxZ60cm",1);
0088 
0089     tree_in -> SetBranchStatus("InttRawHit_bco", 1);
0090     tree_in -> SetBranchStatus("InttRawHit_packetid", 1);
0091     tree_in -> SetBranchStatus("InttRawHit_fee", 1);
0092     tree_in -> SetBranchStatus("InttRawHit_channel_id", 1);
0093     tree_in -> SetBranchStatus("InttRawHit_chip_id", 1);
0094     tree_in -> SetBranchStatus("InttRawHit_adc", 1);
0095     tree_in -> SetBranchStatus("InttRawHit_FPHX_BCO", 1);
0096 
0097     // note : INTT vertex Z
0098     tree_in -> SetBranchStatus("INTTvtxZ", 1);
0099     tree_in -> SetBranchStatus("INTTvtxZError", 1);
0100     tree_in -> SetBranchStatus("NgroupTrapezoidal", 1);
0101     tree_in -> SetBranchStatus("NgroupCoarse", 1);
0102     tree_in -> SetBranchStatus("TrapezoidalFitWidth", 1);
0103     tree_in -> SetBranchStatus("TrapezoidalFWHM", 1);
0104 
0105     tree_in -> SetBranchStatus("NClus", 1);
0106     tree_in -> SetBranchStatus("NClus_Layer1", 1);
0107     tree_in -> SetBranchStatus("ClusAdc", 1);
0108     tree_in -> SetBranchStatus("ClusPhiSize", 1);
0109 
0110 
0111     InttRawHit_bco = 0;
0112     InttRawHit_packetid = 0;
0113     InttRawHit_fee = 0;
0114     InttRawHit_channel_id = 0;
0115     InttRawHit_chip_id = 0;
0116     InttRawHit_adc = 0;
0117     InttRawHit_FPHX_BCO = 0;
0118 
0119     ClusAdc = 0;
0120     ClusPhiSize = 0;
0121 
0122     tree_in -> SetBranchAddress("MBD_z_vtx", &MBD_z_vtx);
0123     tree_in -> SetBranchAddress("is_min_bias", &is_min_bias);
0124     tree_in -> SetBranchAddress("MBD_centrality", &MBD_centrality);
0125     tree_in -> SetBranchAddress("MBD_south_charge_sum", &MBD_south_charge_sum);
0126     tree_in -> SetBranchAddress("MBD_north_charge_sum", &MBD_north_charge_sum);
0127     tree_in -> SetBranchAddress("MBD_charge_sum", &MBD_charge_sum);
0128     tree_in -> SetBranchAddress("MBD_charge_asymm", &MBD_charge_asymm);
0129     tree_in -> SetBranchAddress("InttBcoFullDiff_next", &InttBcoFullDiff_next);
0130 
0131     tree_in -> SetBranchAddress("MBDNSg2",&MBDNSg2);
0132     tree_in -> SetBranchAddress("MBDNSg2_vtxZ10cm",&MBDNSg2_vtxZ10cm);
0133     tree_in -> SetBranchAddress("MBDNSg2_vtxZ30cm",&MBDNSg2_vtxZ30cm);
0134     tree_in -> SetBranchAddress("MBDNSg2_vtxZ60cm",&MBDNSg2_vtxZ60cm);
0135 
0136     tree_in -> SetBranchAddress("InttRawHit_bco", &InttRawHit_bco);
0137     tree_in -> SetBranchAddress("InttRawHit_packetid", &InttRawHit_packetid);
0138     tree_in -> SetBranchAddress("InttRawHit_fee", &InttRawHit_fee);
0139     tree_in -> SetBranchAddress("InttRawHit_channel_id", &InttRawHit_channel_id);
0140     tree_in -> SetBranchAddress("InttRawHit_chip_id", &InttRawHit_chip_id);
0141     tree_in -> SetBranchAddress("InttRawHit_adc", &InttRawHit_adc);
0142     tree_in -> SetBranchAddress("InttRawHit_FPHX_BCO", &InttRawHit_FPHX_BCO);
0143 
0144     // note : INTT vertex Z
0145     tree_in -> SetBranchAddress("INTTvtxZ", &INTTvtxZ);
0146     tree_in -> SetBranchAddress("INTTvtxZError", &INTTvtxZError);
0147     tree_in -> SetBranchAddress("NgroupTrapezoidal", &NgroupTrapezoidal);
0148     tree_in -> SetBranchAddress("NgroupCoarse", &NgroupCoarse);
0149     tree_in -> SetBranchAddress("TrapezoidalFitWidth", &TrapezoidalFitWidth);
0150     tree_in -> SetBranchAddress("TrapezoidalFWHM", &TrapezoidalFWHM);
0151 
0152     tree_in -> SetBranchAddress("NClus", &NClus);
0153     tree_in -> SetBranchAddress("NClus_Layer1", &NClus_Layer1);
0154 
0155     tree_in -> SetBranchAddress("ClusAdc", &ClusAdc);
0156     tree_in -> SetBranchAddress("ClusPhiSize", &ClusPhiSize);
0157 
0158 
0159 }
0160 
0161 void SemiChipClustering::PrepareOutPutFileName()
0162 {
0163     std::string job_index = std::to_string( process_id );
0164     int job_index_len = 5;
0165     job_index.insert(0, job_index_len - job_index.size(), '0');
0166 
0167     std::string runnumber_str = std::to_string( runnumber );
0168     if (runnumber != -1){
0169         int runnumber_str_len = 8;
0170         runnumber_str.insert(0, runnumber_str_len - runnumber_str.size(), '0');
0171     }
0172 
0173     if (output_file_name_suffix.size() > 0 && output_file_name_suffix[0] != '_') {
0174         output_file_name_suffix = "_" + output_file_name_suffix;
0175     }
0176 
0177     output_filename = "SemiChipClus";
0178     output_filename = (runnumber != -1) ? "Data_" + output_filename : "MC_" + output_filename;
0179     
0180     output_filename += (BcoFullDiffCut && runnumber != -1) ? "_BcoFullDiffCut" : "";
0181     output_filename += (INTT_vtxZ_QA) ? "_VtxZQA" : "";
0182 
0183     output_filename += (BadChMask.first) ? "_BadChMask" : "";
0184     output_filename += (ApplyHitQA) ? "_HitQA" : "";
0185     output_filename += (cut_HitBcoDiff.first) ? Form("_HitBcoDiff%d",cut_HitBcoDiff.second) : "";
0186     output_filename += (clone_hit_remove_BCO_tag) ? "_CloneHitRm" : "";
0187 
0188     output_filename += output_file_name_suffix;
0189     output_filename += (runnumber != -1) ? Form("_%s_%s.root",runnumber_str.c_str(),job_index.c_str()) : Form("_%s.root",job_index.c_str());
0190 }
0191 
0192 void SemiChipClustering::PrepareOutPutRootFile()
0193 {
0194     file_out = new TFile(Form("%s/%s", output_directory.c_str(), output_filename.c_str()), "RECREATE");
0195 }
0196 
0197 void SemiChipClustering::PrepareHists()
0198 {
0199     h1D_NSize2Clus = new TH1D("h1D_NSize2Clus","h1D_NSize2Clus;NSize2Clus;Entries", 30, 0, 30);
0200     h1D_NSize1Clus = new TH1D("h1D_NSize1Clus","h1D_NSize1Clus;NSize1Clus;Entries", 45, 0, 45);
0201     
0202     h2D_NClus_NSize2Clus = new TH2D("h2D_NClus_NSize2Clus","h2D_NClus_NSize2Clus;NClus;NSize2Clus", 200, 0, 10000, 30, 0, 30);
0203     h2D_NClus_NSize1Clus = new TH2D("h2D_NClus_NSize1Clus","h2D_NClus_NSize1Clus;NClus;NSize1Clus", 200, 0, 10000, 45, 0, 45);
0204     h2D_nInttRawHit_NSize2Clus = new TH2D("h2D_nInttRawHit_NSize2Clus","h2D_nInttRawHit_NSize2Clus;NInttRawHits;NSize2Clus",200, 0, 25000, 30, 0, 30);
0205     h2D_nInttRawHit_NSize1Clus = new TH2D("h2D_nInttRawHit_NSize1Clus","h2D_nInttRawHit_NSize1Clus;NInttRawHits;NSize1Clus",200, 0, 25000, 45, 0, 45);
0206 
0207     h2D_nChipHit_NSize2Clus = new TH2D("h2D_nChipHit_NSize2Clus","h2D_nChipHit_NSize2Clus;NChipHits;NSize2Clus",128, 0, 128, 30, 0, 30);
0208     h2D_nChipHit_NSize1Clus = new TH2D("h2D_nChipHit_NSize1Clus","h2D_nChipHit_NSize1Clus;NChipHits;NSize1Clus",128, 0, 128, 45, 0, 45);
0209 
0210     h2D_NSize2Clus_NSize1Clus = new TH2D("h2D_NSize2Clus_NSize1Clus","h2D_NSize2Clus_NSize1Clus;NSize2Clus;NSize1Clus", 30, 0, 30, 45, 0, 45);
0211     h2D_NSize2Clus_LargestSize = new TH2D("h2D_NSize2Clus_LargestSize","h2D_NSize2Clus_LargestSize;NSize2Clus;LargestSize", 30, 0, 30, 128, 0, 128);
0212 
0213     h1D_ClusSizeCount = new TH1D("h1D_ClusSizeCount","h1D_ClusSizeCount;ClusSize;Entries", 130,0,130);
0214     h1D_ClusSizeCount_all = new TH1D("h1D_ClusSizeCount_all","h1D_ClusSizeCount_all;ClusSize;Entries", 130,0,130);
0215 
0216     h1D_confirm_HitBcoDiff_post = new TH1D("h1D_confirm_HitBcoDiff_post","h1D_confirm_HitBcoDiff_post;Time_bucket;Entries",128, 0, 128);
0217     h1D_confirm_HitBcoDiff = new TH1D("h1D_confirm_HitBcoDiff","h1D_confirm_HitBcoDiff;Time_bucket;Entries",128, 0, 128);
0218     h2D_confirm_NClus_MBDChargeSum = new TH2D("h2D_confirm_NClus_MBDChargeSum","h2D_confirm_NClus_MBDChargeSum;INTT NClus; MBD Charge Sum",200, 0, 10000, 200, 0, 3000);
0219 
0220     h2D_nChipHit_Size2DistAvg = new TH2D("h2D_nChipHit_Size2DistAvg","h2D_nChipHit_Size2DistAvg;NChipHits;Avg. Size2Dist",128, 0, 128, 100, 0, 20);
0221     h2D_nChipHit_Size2DistStd = new TH2D("h2D_nChipHit_Size2DistStd","h2D_nChipHit_Size2DistStd;NChipHits;StdDev Size2Dist",128, 0, 128, 200, 0, 20);
0222     h2D_NSize2Clus_Size2Dist = new TH2D("h2D_NSize2Clus_Size2Dist","h2D_NSize2Clus_Size2Dist;NSize2Clus;Avg. Size2Dist",30,0,30,100,0,20);
0223     h2D_Size2DistAvg_Size2DistStd = new TH2D("h2D_Size2DistAvg_Size2DistStd","h2D_Size2DistAvg_Size2DistStd;Avg. Size2Dist;StdDev Size2Dist",100,0,20,200,0,20);
0224     h2D_Size2DistAvg_Size2DistStd_WithLargeSize = new TH2D("h2D_Size2DistAvg_Size2DistStd_WithLargeSize","h2D_Size2DistAvg_Size2DistStd_WithLargeSize;Avg. Size2Dist;StdDev Size2Dist",100,0,20,200,0,20);
0225 
0226 
0227     h2D_Inclusive_Inner_Outer_NClus = new TH2D("h2D_Inclusive_Inner_Outer_NClus","h2D_Inclusive_Inner_Outer_NClus;NClus (Inner);NClus (Outer)",200,0,5000,200,0,5000);
0228     h2D_Inclusive_NClus_MBDChargeSum = new TH2D("h2D_Inclusive_NClus_MBDChargeSum","h2D_Inclusive_NClus_MBDChargeSum;NClus;MBD charge sum",200,0,10000,200,0,3000);
0229     h1D_Inclusive_MBDCentrality = new TH1D("h1D_Inclusive_MBDCentrality","h1D_Inclusive_MBDCentrality;MBD Centrality [%];Entries",nCentralityFineBin, CentralityFineEdge_min, CentralityFineEdge_max);
0230     h1D_Inclusive_ClusPhiSize = new TH1D("h1D_Inclusive_ClusPhiSize","h1D_Inclusive_ClusPhiSize;ClusPhiSize;Entries",128,0,128);
0231     h1D_Inclusive_ClusAdc = new TH1D("h1D_Inclusive_ClusAdc","h1D_Inclusive_ClusAdc;ClusAdc;Entries",200,0,18000);
0232 
0233     h2D_Post_Inner_Outer_NClus = new TH2D("h2D_Post_Inner_Outer_NClus","h2D_Post_Inner_Outer_NClus;NClus (Inner);NClus (Outer)",200,0,5000,200,0,5000);
0234     h2D_Post_NClus_MBDChargeSum = new TH2D("h2D_Post_NClus_MBDChargeSum","h2D_Post_NClus_MBDChargeSum;NClus;MBD charge sum",200,0,10000,200,0,3000);
0235     h1D_Post_MBDCentrality = new TH1D("h1D_Post_MBDCentrality","h1D_Post_MBDCentrality;MBD Centrality [%];Entries",nCentralityFineBin, CentralityFineEdge_min, CentralityFineEdge_max);
0236     h1D_Post_ClusPhiSize = new TH1D("h1D_Post_ClusPhiSize","h1D_Post_ClusPhiSize;ClusPhiSize;Entries",128,0,128);
0237     h1D_Post_ClusAdc = new TH1D("h1D_Post_ClusAdc","h1D_Post_ClusAdc;ClusAdc;Entries",200,0,18000);
0238 
0239     h2D_Killed_Inner_Outer_NClus = new TH2D("h2D_Killed_Inner_Outer_NClus","h2D_Killed_Inner_Outer_NClus;NClus (Inner);NClus (Outer)",200,0,5000,200,0,5000);
0240     h2D_Killed_NClus_MBDChargeSum = new TH2D("h2D_Killed_NClus_MBDChargeSum","h2D_Killed_NClus_MBDChargeSum;NClus;MBD charge sum",200,0,10000,200,0,3000);
0241     h1D_Killed_MBDCentrality = new TH1D("h1D_Killed_MBDCentrality","h1D_Killed_MBDCentrality;MBD Centrality [%];Entries",nCentralityFineBin, CentralityFineEdge_min, CentralityFineEdge_max);
0242     h1D_Killed_ClusPhiSize = new TH1D("h1D_Killed_ClusPhiSize","h1D_Killed_ClusPhiSize;ClusPhiSize;Entries",128,0,128);
0243     h1D_Killed_ClusAdc = new TH1D("h1D_Killed_ClusAdc","h1D_Killed_ClusAdc;ClusAdc;Entries",200,0,18000);
0244 
0245 
0246     h1D_chip_hit_container_map.clear();
0247     for (int server_i = 0; server_i < nFelix; server_i++){
0248         for (int Fch_i = 0; Fch_i < nFelix_channel; Fch_i++){
0249             for (int chip_i = 0; chip_i < nChip; chip_i++){
0250                 
0251                 h1D_chip_hit_container_map[Form("server%d_Fch%d_chip%d",server_i,Fch_i,chip_i)] = 
0252                 new TH1D(Form("server%d_Fch%d_chip%d",server_i,Fch_i,chip_i), Form("server%d_Fch%d_chip%d",server_i,Fch_i,chip_i), 128, 0, 128);
0253 
0254             }
0255         }
0256     }
0257 }
0258 
0259 void SemiChipClustering::EvtCleanUp()
0260 {
0261     evt_inttHits_map.clear();
0262 
0263     for ( auto &pair : h1D_chip_hit_container_map)
0264     {
0265         pair.second -> Reset("ICESM");
0266     }
0267 
0268     h1D_ClusSizeCount -> Reset("ICESM");
0269 }
0270 
0271 bool SemiChipClustering::DoSemiChipCluster(
0272     int eID_count,
0273     int NClus_in,
0274     std::vector<unsigned long> bco_vec_in,
0275     std::vector<unsigned int> packetid_vec_in,
0276     std::vector<unsigned short> fee_vec_in,
0277     std::vector<unsigned short> channel_id_vec_in,
0278     std::vector<unsigned short> chip_id_vec_in,
0279     std::vector<unsigned short> adc_vec_in,
0280     std::vector<unsigned short> FPHX_BCO_vec_in
0281 )
0282 {
0283     EvtCleanUp();
0284 
0285     for (int hit_i = 0; hit_i < bco_vec_in.size(); hit_i++)
0286     {
0287         if (hit_i == 0) {evt_INTT_bco_full = bco_vec_in[hit_i];}
0288 
0289         int bco_full = bco_vec_in[hit_i];
0290 
0291         int server = packetid_vec_in[hit_i] - Felix_offset; // note : the felix server ID
0292         int felix_ch = fee_vec_in[hit_i]; // note : the felix channel ID 0 - 13
0293         int chip = (chip_id_vec_in[hit_i] - 1) % 26; // note : chip ID 0 - 25
0294         int channel = channel_id_vec_in[hit_i]; // note : channel ID 0 - 127
0295         int adc = adc_vec_in[hit_i]; // note : adc value
0296         int bco = FPHX_BCO_vec_in[hit_i];
0297         int bco_diff = (bco - (bco_full & 0x7fU) + 128) % 128;
0298 
0299         h1D_confirm_HitBcoDiff -> Fill(bco_diff);
0300 
0301 
0302         if (eID_count % 1000 == 0 && hit_i%25 == 0) {
0303             std::cout << "eID: "<< eID_count <<" INTTBcoResolution::PrepareINTT - server: "<< server << " felix_ch: "<< felix_ch << " chip: "<< chip << " channel: "<< channel << " adc: "<< adc << " bco: "<< bco << " bco_diff: "<< bco_diff << std::endl;
0304         }
0305 
0306         if (ApplyHitQA && (server < 0 || server > 7)) // note : server ID 0 - 7
0307         {
0308             std::cout << "eID: "<< eID_count <<" INTTBcoResolution::PrepareINTT - server out of range, this server is : "<< server << std::endl;
0309             continue;
0310         }
0311 
0312         if (ApplyHitQA && (felix_ch < 0 || felix_ch > 13)) // note : felix channel ID 0 - 13
0313         {
0314             std::cout << "eID: "<< eID_count <<" INTTBcoResolution::PrepareINTT - felix channel out of range, this channel is : "<< felix_ch << std::endl;
0315             continue;
0316         }
0317 
0318         if (ApplyHitQA && (chip < 0 || chip > 25)) // note : chip ID 0 - 25
0319         {
0320             std::cout << "eID: "<< eID_count <<" INTTBcoResolution::PrepareINTT - chip out of range, this chip is : "<< chip << std::endl;
0321             continue;
0322         }
0323 
0324         if (ApplyHitQA && (channel < 0 || channel > 127)) // note : channel ID 0 - 127
0325         {
0326             std::cout << "eID: "<< eID_count <<" INTTBcoResolution::PrepareINTT - channel out of range, this channel is : "<< channel << std::endl;
0327             continue;
0328         }
0329 
0330         if (ApplyHitQA && (adc < 0 || adc > 7)) // note : adc ID 0 - 7
0331         {
0332             std::cout << "eID: "<< eID_count <<" INTTBcoResolution::PrepareINTT - adc out of range, this adc is : "<< adc << std::endl;
0333             continue;
0334         }
0335         
0336         if (ApplyHitQA && (bco < 0 || bco > 127)) // note : bco ID 0 - 127
0337         {
0338             std::cout << "eID: "<< eID_count <<" INTTBcoResolution::PrepareINTT - bco out of range, this bco is : "<< bco << std::endl;
0339             continue;
0340         }
0341 
0342         if (BadChMask.first && (hot_channel_map.find(Form("%d_%d_%d_%d",server,felix_ch,chip,channel)) != hot_channel_map.end())) // note : if we want to remove the hot channels
0343         {
0344             continue;
0345         }
0346 
0347         if (cut_HitBcoDiff.first && (bco_diff < cut_HitBcoDiff.second - 1 || bco_diff > cut_HitBcoDiff.second + 1)) // note : bco_diff cut
0348         {
0349             // std::cout << "eID: "<< eID_count <<" INTTBcoResolution::PrepareINTT - bco_diff out of range, this bco_diff is : "<< bco_diff << std::endl;
0350             continue;
0351         }
0352 
0353         SemiChipClustering::inttHitstr this_hit;
0354         
0355         this_hit.hit_server = server;
0356         this_hit.hit_felix_ch = felix_ch;
0357         this_hit.hit_chip = chip;
0358         this_hit.hit_channel = channel;
0359         this_hit.hit_adc = adc_conversion_vec[adc];
0360         this_hit.hit_bco = bco;
0361         this_hit.hit_bco_diff = bco_diff;
0362 
0363         if (clone_hit_remove_BCO_tag) // note : can only be one of each
0364         {
0365             if (evt_inttHits_map.find(Form("%d_%d_%d_%d_%d",server,felix_ch,chip,channel,bco)) == evt_inttHits_map.end()) // note : if it's not found, we just add it
0366             {
0367                 evt_inttHits_map[Form("%d_%d_%d_%d_%d",server,felix_ch,chip,channel,bco)] = this_hit;
0368             }
0369         }
0370         else // note : if we don't want to remove the clone hits
0371         {
0372             evt_inttHits_map[Form("%d",hit_i)] = this_hit; // note : only index i to make the key unique
0373         }
0374 
0375     } // note : end of loop over hits
0376 
0377     for (auto pair : evt_inttHits_map)
0378     {
0379         SemiChipClustering::inttHitstr this_hit = pair.second;
0380         int server = this_hit.hit_server;
0381         int felix_ch = this_hit.hit_felix_ch;
0382         int chip = this_hit.hit_chip;
0383         int channel = this_hit.hit_channel;
0384         int adc = this_hit.hit_adc;
0385         int bco = this_hit.hit_bco;
0386         int bco_diff = this_hit.hit_bco_diff;
0387 
0388         h1D_confirm_HitBcoDiff_post -> Fill(bco_diff);
0389 
0390         h1D_chip_hit_container_map[Form("server%d_Fch%d_chip%d",server,felix_ch,chip)] -> Fill(channel, adc);
0391     }
0392 
0393     bool good_evt_flag = true;
0394 
0395     for (auto pair : h1D_chip_hit_container_map)
0396     {
0397         if (pair.second -> GetEntries() == 0) {continue;}
0398         
0399         SemiChipClustering::chip_clu_info clu_info = find_Ngroup(pair.second);
0400         h1D_ClusSizeCount -> Reset("ICESM");
0401 
0402         for (auto size : clu_info.size_vec)
0403         {
0404             h1D_ClusSizeCount -> Fill(size);
0405             h1D_ClusSizeCount_all -> Fill(size);
0406         }
0407 
0408         int NSize1Clus = h1D_ClusSizeCount -> GetBinContent(2);
0409         int NSize2Clus = h1D_ClusSizeCount -> GetBinContent(3);
0410 
0411         h1D_NSize1Clus -> Fill(NSize1Clus);
0412         h2D_NClus_NSize1Clus -> Fill(NClus_in, NSize1Clus);
0413         h2D_nInttRawHit_NSize1Clus -> Fill(evt_inttHits_map.size(), NSize1Clus);
0414         h2D_nChipHit_NSize1Clus -> Fill(pair.second -> GetEntries(), NSize1Clus);
0415 
0416         h2D_NSize2Clus_NSize1Clus -> Fill(NSize2Clus, NSize1Clus);
0417         h2D_NSize2Clus_LargestSize -> Fill(NSize2Clus, clu_info.largest_size);
0418 
0419 
0420         std::pair<double,double> size2dist_pair = GetSize2Dist(clu_info, eID_count, pair.first);
0421         if (size2dist_pair.first == size2dist_pair.first && size2dist_pair.second == size2dist_pair.second){
0422             double avg_size2dist = size2dist_pair.first;
0423             double std_size2dist = size2dist_pair.second;
0424 
0425             h2D_nChipHit_Size2DistAvg -> Fill(pair.second -> GetEntries(), avg_size2dist);
0426             h2D_nChipHit_Size2DistStd -> Fill(pair.second -> GetEntries(), std_size2dist);
0427 
0428             h2D_NSize2Clus_Size2Dist -> Fill(NSize2Clus, avg_size2dist);
0429             h2D_Size2DistAvg_Size2DistStd -> Fill(avg_size2dist, std_size2dist);
0430 
0431             // if (avg_size2dist >= 4 && avg_size2dist < 5 && std_size2dist < 0.3 && pair.second -> GetEntries() > 68){ // todo : the chip hit multiplicity cut
0432             //     good_evt_flag = false;
0433             // } 
0434 
0435             // if (pair.second -> GetEntries() > 60){
0436             //     good_evt_flag = false;
0437             // }
0438 
0439             if (clu_info.largest_size > 20 && NSize2Clus > 10){
0440                 good_evt_flag = false;
0441             }
0442             
0443             if (clu_info.largest_size <= 20) {continue;}    
0444 
0445             h2D_Size2DistAvg_Size2DistStd_WithLargeSize -> Fill(avg_size2dist, std_size2dist);
0446         }
0447 
0448         if (clu_info.largest_size <= 20) {continue;}
0449 
0450         h1D_NSize2Clus -> Fill(NSize2Clus);        
0451         h2D_NClus_NSize2Clus -> Fill(NClus_in, NSize2Clus);        
0452         h2D_nInttRawHit_NSize2Clus -> Fill(evt_inttHits_map.size(), NSize2Clus);        
0453         h2D_nChipHit_NSize2Clus -> Fill(pair.second -> GetEntries(), NSize2Clus);    
0454     }
0455 
0456     return good_evt_flag;
0457 }
0458 
0459 std::pair<double, double> SemiChipClustering::GetSize2Dist(SemiChipClustering::chip_clu_info clu_info_in, int eID_count, std::string chip_string)
0460 {
0461     std::vector<std::pair<double, double>> obj_vec; obj_vec.clear();
0462 
0463     for (int clu_i = 0; clu_i < clu_info_in.size_vec.size(); clu_i++){
0464         if (clu_info_in.size_vec[clu_i] != 2) {continue;}
0465 
0466         obj_vec.push_back( std::make_pair(clu_info_in.edge_l_vec[clu_i], clu_info_in.edge_r_vec[clu_i]) );
0467     }
0468 
0469     // note : no cluster with size 2 or only 1 cluster with size 2
0470     if (obj_vec.size() < 2) {return std::make_pair(std::nan(""),std::nan(""));}
0471 
0472     if (eID_count % 200 == 0){
0473         for (int i = 0; i < obj_vec.size(); i++){
0474             std::cout<<"In GetSize2Dist, eID: "<<eID_count<<", "<<chip_string<<", edge : "<<obj_vec[i].first<<" "<<obj_vec[i].second<<std::endl;
0475         }
0476     }
0477 
0478     std::vector<double> distances; distances.clear();
0479 
0480     for (size_t i = 0; i < obj_vec.size() - 1; ++i) {
0481         double distance = obj_vec[i + 1].first - obj_vec[i].second;
0482         distances.push_back(distance);
0483     }
0484 
0485     return std::make_pair(vector_average(distances), vector_stddev(distances));
0486 
0487 }
0488 
0489 double  SemiChipClustering::vector_average (std::vector <double> input_vector) {
0490     return accumulate( input_vector.begin(), input_vector.end(), 0.0 ) / double(input_vector.size());
0491 }
0492 
0493 double SemiChipClustering::vector_stddev (std::vector <double> input_vector){
0494     
0495     double sum_subt = 0;
0496     double average  = accumulate( input_vector.begin(), input_vector.end(), 0.0 ) / double(input_vector.size());
0497     
0498     // cout<<"average is : "<<average<<endl;
0499 
0500     for (int i=0; i<int(input_vector.size()); i++){ sum_subt += pow((input_vector[i] - average),2); }
0501 
0502     //cout<<"sum_subt : "<<sum_subt<<endl;
0503     // cout<<"print from the function, average : "<<average<<" std : "<<stddev<<endl;
0504 
0505     return sqrt( sum_subt / double(input_vector.size()-1) );
0506 }   
0507 
0508 SemiChipClustering::chip_clu_info SemiChipClustering::find_Ngroup(TH1D * hist_in)
0509 {
0510     double Highest_bin_Content __attribute__((unused)) = hist_in -> GetBinContent(hist_in -> GetMaximumBin());
0511     double Highest_bin_Center   = hist_in -> GetBinCenter(hist_in -> GetMaximumBin());
0512 
0513     int group_Nbin = 0;
0514     int peak_group_ID = -9999;
0515     double group_entry = 0;
0516     double peak_group_ratio __attribute__((unused));
0517     std::vector<int> group_Nbin_vec; group_Nbin_vec.clear();
0518     std::vector<double> group_entry_vec; group_entry_vec.clear();
0519     std::vector<double> group_widthL_vec; group_widthL_vec.clear();
0520     std::vector<double> group_widthR_vec; group_widthR_vec.clear();
0521     std::vector<double> clu_pos_vec; clu_pos_vec.clear();
0522 
0523     for (int i = 0; i < hist_in -> GetNbinsX(); i++){
0524         // todo : the background rejection is here : Highest_bin_Content/2. for the time being
0525         // double bin_content = ( hist_in -> GetBinContent(i+1) <= Highest_bin_Content/2.) ? 0. : ( hist_in -> GetBinContent(i+1) - Highest_bin_Content/2. );
0526 
0527         double bin_content = hist_in -> GetBinContent(i+1);
0528 
0529         if (bin_content != 0){
0530             
0531             if (group_Nbin == 0) {
0532                 group_widthL_vec.push_back(hist_in -> GetBinCenter(i+1) - (hist_in -> GetBinWidth(i+1)/2.));
0533             }
0534 
0535             group_Nbin += 1; 
0536             group_entry += bin_content;
0537         }
0538         else if (bin_content == 0 && group_Nbin != 0){
0539             group_widthR_vec.push_back(hist_in -> GetBinCenter(i+1) - (hist_in -> GetBinWidth(i+1)/2.));
0540             group_Nbin_vec.push_back(group_Nbin);
0541             group_entry_vec.push_back(group_entry);
0542             group_Nbin = 0;
0543             group_entry = 0;
0544         }
0545     }
0546     if (group_Nbin != 0) {
0547         group_Nbin_vec.push_back(group_Nbin);
0548         group_entry_vec.push_back(group_entry);
0549         group_widthR_vec.push_back(hist_in -> GetXaxis()->GetXmax());
0550     } // note : the last group at the edge
0551 
0552     // note : find the peak group
0553     for (int i = 0; i < int(group_Nbin_vec.size()); i++){
0554         if (group_widthL_vec[i] < Highest_bin_Center && Highest_bin_Center < group_widthR_vec[i]){
0555             peak_group_ID = i;
0556             break;
0557         }
0558     }
0559     
0560     // note : On Nov 6, 2024, we no longer need to calculate the ratio of the peak group
0561     // double denominator = (accumulate( group_entry_vec.begin(), group_entry_vec.end(), 0.0 ));
0562     // denominator = (denominator <= 0) ? 1. : denominator;
0563     // peak_group_ratio = group_entry_vec[peak_group_ID] / denominator;
0564     peak_group_ratio = -999;
0565 
0566     double peak_group_left __attribute__((unused))  = (double(group_Nbin_vec.size()) == 0) ? -999 : group_widthL_vec[peak_group_ID];
0567     double peak_group_right __attribute__((unused)) = (double(group_Nbin_vec.size()) == 0) ? 999  : group_widthR_vec[peak_group_ID];
0568 
0569     for (int i = 0; i < int(group_Nbin_vec.size()); i++){
0570         clu_pos_vec.push_back( CoM_cluster_pos(hist_in, group_widthL_vec[i], group_widthR_vec[i]) );
0571     }
0572 
0573     SemiChipClustering::chip_clu_info out_str;
0574 
0575     out_str.adc_vec = group_entry_vec;
0576     out_str.size_vec = group_Nbin_vec;
0577     out_str.edge_l_vec = group_widthL_vec;
0578     out_str.edge_r_vec = group_widthR_vec;
0579     out_str.pos_vec = clu_pos_vec;
0580     out_str.largest_size = *(std::max_element(group_Nbin_vec.begin(), group_Nbin_vec.end()));
0581 
0582     // for (int i = 0; i < group_Nbin_vec.size(); i++)
0583     // {
0584     //     std::cout<<" "<<std::endl;
0585     //     std::cout<<"group width : "<<group_Nbin_vec[i]<<std::endl;
0586     //     std::cout<<"group adc : "<<group_entry_vec[i]<<std::endl;
0587     //     std::cout<<group_widthL_vec[i]<<" "<<group_widthR_vec[i]<<std::endl;
0588     // }
0589 
0590     // cout<<" "<<endl;
0591     // cout<<"N group : "<<group_Nbin_vec.size()<<endl;
0592     // cout<<"Peak group ID : "<<peak_group_ID<<endl;
0593     // cout<<"peak group width : "<<group_widthL_vec[peak_group_ID]<<" "<<group_widthR_vec[peak_group_ID]<<endl;
0594     // cout<<"ratio : "<<peak_group_ratio<<endl;
0595     
0596     // note : {N_group, ratio (if two), peak widthL, peak widthR}
0597     return out_str;
0598 }
0599 
0600 double SemiChipClustering::CoM_cluster_pos(TH1D * hist_in, double edge_l, double edge_r)
0601 {
0602     double sum = 0;
0603     double sum_weight = 0;
0604 
0605     for (int i = 0; i < hist_in -> GetNbinsX(); i++){
0606         double bin_content = hist_in -> GetBinContent(i+1);
0607         double bin_center = hist_in -> GetBinCenter(i+1);
0608 
0609         if (edge_l < bin_center && bin_center < edge_r){
0610             sum += bin_content * bin_center;
0611             sum_weight += bin_content;
0612         }
0613     }
0614 
0615     return (sum_weight == 0) ? -999 : (sum / sum_weight) - hist_in -> GetBinWidth(1)/2.;
0616 }
0617 
0618 void SemiChipClustering::MainProcess()
0619 {
0620     for (int i = 0; i < run_nEvents; i++)
0621     {
0622         tree_in -> GetEntry(i);
0623         EvtCleanUp();
0624 
0625         if (i % 10 == 0) {std::cout << "Processing event " << i<<", NClus : "<< NClus << std::endl;}
0626 
0627         // =======================================================================================================================================================
0628         // note : hard cut
0629 
0630         // note : for data
0631         if (runnumber != -1 && BcoFullDiffCut && InttBcoFullDiff_next <= cut_InttBcoFullDIff_next) {continue;}
0632         if (runnumber != -1 && MBDNSg2 != 1) {continue;} // todo: assume MC no trigger
0633 
0634         // note : for MC
0635         // if (runnumber == -1 && NTruthVtx != 1) {continue;}
0636 
0637         // note : both data and MC
0638         if (is_min_bias != 1) {continue;}
0639         if (MBD_z_vtx != MBD_z_vtx) {continue;}
0640         if (MBD_centrality != MBD_centrality) {continue;}
0641         if (MBD_centrality < 0 || MBD_centrality > 1) {continue;}
0642         if (INTTvtxZ != INTTvtxZ) {continue;}
0643         if (MBD_z_vtx < cut_GlobalMBDvtxZ.first || MBD_z_vtx > cut_GlobalMBDvtxZ.second) {continue;} // todo: the hard cut 60 cm 
0644         // =======================================================================================================================================================
0645 
0646         // =======================================================================================================================================================
0647         // note : optional cut
0648         if (INTT_vtxZ_QA && (MBD_z_vtx - INTTvtxZ < cut_vtxZDiff.first || MBD_z_vtx - INTTvtxZ > cut_vtxZDiff.second) ) {continue;}
0649         if (INTT_vtxZ_QA && (TrapezoidalFitWidth < cut_TrapezoidalFitWidth.first || TrapezoidalFitWidth > cut_TrapezoidalFitWidth.second)){continue;}
0650         if (INTT_vtxZ_QA && (TrapezoidalFWHM < cut_TrapezoidalFWHM.first || TrapezoidalFWHM > cut_TrapezoidalFWHM.second)){continue;}
0651         if (INTT_vtxZ_QA && (INTTvtxZError < cut_INTTvtxZError.first || INTTvtxZError > cut_INTTvtxZError.second)){continue;}
0652         // =======================================================================================================================================================
0653 
0654         h2D_confirm_NClus_MBDChargeSum -> Fill(NClus, MBD_charge_sum);
0655 
0656         bool good_evt_flag = DoSemiChipCluster(
0657             i,
0658             NClus,
0659             *InttRawHit_bco,
0660             *InttRawHit_packetid,
0661             *InttRawHit_fee,
0662             *InttRawHit_channel_id,
0663             *InttRawHit_chip_id,
0664             *InttRawHit_adc,
0665             *InttRawHit_FPHX_BCO
0666         );
0667 
0668         TH2D * temp_h2D_Inner_Outer_NClus = (good_evt_flag) ? h2D_Post_Inner_Outer_NClus : h2D_Killed_Inner_Outer_NClus;
0669         TH2D * temp_h2D_NClus_MBDChargeSum = (good_evt_flag) ? h2D_Post_NClus_MBDChargeSum : h2D_Killed_NClus_MBDChargeSum;
0670         TH1D * temp_h1D_MBDCentrality = (good_evt_flag) ? h1D_Post_MBDCentrality : h1D_Killed_MBDCentrality;
0671         TH1D * temp_h1D_ClusPhiSize = (good_evt_flag) ? h1D_Post_ClusPhiSize : h1D_Killed_ClusPhiSize;
0672         TH1D * temp_h1D_ClusAdc = (good_evt_flag) ? h1D_Post_ClusAdc : h1D_Killed_ClusAdc;
0673         
0674         h2D_Inclusive_Inner_Outer_NClus -> Fill(NClus_Layer1,NClus - NClus_Layer1);
0675         h2D_Inclusive_NClus_MBDChargeSum -> Fill(NClus,MBD_charge_sum);
0676         h1D_Inclusive_MBDCentrality -> Fill(MBD_centrality);
0677 
0678         temp_h2D_Inner_Outer_NClus -> Fill(NClus_Layer1,NClus - NClus_Layer1);
0679         temp_h2D_NClus_MBDChargeSum -> Fill(NClus,MBD_charge_sum);
0680         temp_h1D_MBDCentrality -> Fill(MBD_centrality);
0681 
0682         for (int clu_i = 0; clu_i < ClusPhiSize -> size(); clu_i++)
0683         {
0684             h1D_Inclusive_ClusPhiSize -> Fill(ClusPhiSize -> at(clu_i));
0685             h1D_Inclusive_ClusAdc -> Fill(ClusAdc -> at(clu_i));
0686 
0687             temp_h1D_ClusPhiSize -> Fill(ClusPhiSize -> at(clu_i));
0688             temp_h1D_ClusAdc -> Fill(ClusAdc -> at(clu_i));
0689         }
0690     }
0691 }
0692 
0693 void SemiChipClustering::EndRun()
0694 {
0695     file_out -> cd();
0696 
0697     h1D_confirm_HitBcoDiff -> Write();
0698     h1D_confirm_HitBcoDiff_post -> Write();
0699     h2D_confirm_NClus_MBDChargeSum -> Write();
0700 
0701     h1D_ClusSizeCount_all -> Write();
0702     h1D_NSize2Clus -> Write();
0703     h1D_NSize1Clus -> Write();
0704     h2D_NClus_NSize2Clus -> Write();
0705     h2D_NClus_NSize1Clus -> Write();
0706     h2D_nInttRawHit_NSize2Clus -> Write();
0707     h2D_nInttRawHit_NSize1Clus -> Write();
0708     h2D_nChipHit_NSize2Clus -> Write();
0709     h2D_nChipHit_NSize1Clus -> Write();
0710 
0711     h2D_NSize2Clus_NSize1Clus -> Write();
0712     h2D_NSize2Clus_LargestSize -> Write();
0713 
0714     h2D_nChipHit_Size2DistStd -> Write();
0715     h2D_nChipHit_Size2DistAvg -> Write();
0716     h2D_NSize2Clus_Size2Dist -> Write();
0717     h2D_Size2DistAvg_Size2DistStd -> Write();
0718     h2D_Size2DistAvg_Size2DistStd_WithLargeSize -> Write();
0719 
0720     h2D_Inclusive_Inner_Outer_NClus -> Write();
0721     h2D_Inclusive_NClus_MBDChargeSum -> Write();
0722     h1D_Inclusive_MBDCentrality -> Write();
0723     h1D_Inclusive_ClusPhiSize -> Write();
0724     h1D_Inclusive_ClusAdc -> Write();
0725     
0726     h2D_Post_Inner_Outer_NClus -> Write();
0727     h2D_Post_NClus_MBDChargeSum -> Write();
0728     h1D_Post_MBDCentrality -> Write();
0729     h1D_Post_ClusPhiSize -> Write();
0730     h1D_Post_ClusAdc -> Write();
0731     
0732     h2D_Killed_Inner_Outer_NClus -> Write();
0733     h2D_Killed_NClus_MBDChargeSum -> Write();
0734     h1D_Killed_MBDCentrality -> Write();
0735     h1D_Killed_ClusPhiSize -> Write();
0736     h1D_Killed_ClusAdc -> Write();
0737 
0738     file_out -> Close();
0739 }
0740 
0741 void SemiChipClustering::PrepareHotChannel()
0742 {
0743   hot_file_in = TFile::Open(BadChMask.second.c_str());
0744 
0745   if (!hot_file_in)
0746   {
0747     std::cout << "INTTBcoResolution::PrepareHotChannel - hot channel file not found" << std::endl;
0748     exit(1);
0749   }
0750 
0751   hot_tree_in = (TTree*)hot_file_in->Get(hot_tree_name.c_str());
0752 
0753   hot_tree_in->SetBranchAddress("IID",&IID);
0754   hot_tree_in->SetBranchAddress("Ichannel",&Ichannel);
0755   hot_tree_in->SetBranchAddress("Ichip",&Ichip);
0756   hot_tree_in->SetBranchAddress("Ifelix_channel",&Ifelix_channel);
0757   hot_tree_in->SetBranchAddress("Ifelix_server",&Ifelix_server);
0758   hot_tree_in->SetBranchAddress("Iflag",&Iflag);
0759 
0760   hot_channel_map.clear();
0761 
0762   for (int i=0; i<hot_tree_in->GetEntries(); i++)
0763   {
0764     hot_tree_in->GetEntry(i);
0765     hot_channel_map[Form("%d_%d_%d_%d",Ifelix_server,Ifelix_channel,Ichip,Ichannel)] = Form("%d",Iflag);
0766   }
0767 
0768   return;
0769 }
[ 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!