mvtx/MvtxHotPixel/MvtxFakeHitRate.cc

0001 #include "MvtxFakeHitRate.h"
0002
0003 #include <mvtx/MvtxHitMap.h>
0004 #include <mvtx/MvtxPixelDefs.h>
0005 #include <mvtx/MvtxPixelMask.h>
0006
0007 #include <ffarawobjects/MvtxRawEvtHeader.h>
0008 #include <ffarawobjects/MvtxRawHit.h>
0009 #include <ffarawobjects/MvtxRawHitContainer.h>
0010
0011 #include <fun4all/Fun4AllReturnCodes.h>
0012 #include <fun4all/PHTFileServer.h>
0013
0014 #include <phool/PHCompositeNode.h>
0015 #include <phool/getClass.h>
0016 #include <phool/phool.h>
0017
0018 #include <TH1.h>
0019 #include <TTree.h>
0020
0021 #include <algorithm>
0022 #include <cstdint>
0023 #include <cstdlib>
0024 #include <iostream>
0025 #include <string>
0026 #include <utility>
0027 #include <vector>
0028
0029 // MvtxFakeHitRate class
0030 //==============================================================================
0031 MvtxFakeHitRate::~MvtxFakeHitRate()
0032 {
0033   // clean up
0034
0035
0036     delete m_hot_pixel_mask;
0037
0038
0039
0040     delete m_hit_map;
0041
0042 }
0043
0044 int MvtxFakeHitRate::InitRun(PHCompositeNode* /*topNode*/)
0045 {
0046   // Load the hot pixel map from the CDB
0047   delete m_hot_pixel_mask;  // make cppcheck happy
0048   m_hot_pixel_mask = new MvtxPixelMask();
0049   if (m_load_from_cdb)
0050   {
0051     std::cout << "MvtxFakeHitRate::InitRun - Loading hot pixel map from CDB" << std::endl;
0052     m_hot_pixel_mask->load_from_CDB();
0053   }
0054
0055   // Create the hit map
0056   m_hit_map = new MvtxHitMap();
0057   m_hit_map->clear();
0058
0059   std::cout << "MvtxFakeHitRate::InitRun - Writing output to " << m_outputfile << std::endl;
0060
0061   // Create the output file and trees
0062   PHTFileServer::open(m_outputfile, "RECREATE");
0063   // main tree
0064   m_tree = new TTree("masked_pixels", "masked_pixels");
0065   m_tree->OptimizeBaskets();
0066   m_tree->SetAutoSave(-5e6);
0067   m_tree->Branch("num_strobes", &m_num_strobes);
0068   m_tree->Branch("num_masked_pixels", &m_num_masked_pixels);
0069   m_tree->Branch("noise_threshold", &m_noise_threshold);
0070   m_tree->Branch("masked_pixels", &m_masked_pixels);
0071
0072   // current mask tree
0073   m_current_mask = new TTree("current_mask", "current_mask");
0074   m_current_mask->OptimizeBaskets();
0075   m_current_mask->SetAutoSave(-5e6);
0076   m_current_mask->Branch("nmasked", &m_current_nmasked);
0077   m_current_mask->Branch("threshold", &m_current_threshold);
0078   m_current_mask->Branch("masked_pixels", &m_current_masked_pixels);
0079
0080   // threshold vs nmasked histogram
0081   m_threshold_vs_nmasked = new TH1D("threshold_vs_nmasked", "threshold_vs_nmasked", m_max_masked_pixels, 0.5, (1.0 * m_max_masked_pixels) + 0.5);
0082
0083   m_num_strobes = 0;
0084   return Fun4AllReturnCodes::EVENT_OK;
0085 }
0086
0087 int MvtxFakeHitRate::get_nodes(PHCompositeNode* topNode)
0088 {
0089   // get dst nodes
0090   m_mvtx_raw_event_header = findNode::getClass<MvtxRawEvtHeader>(topNode, "MVTXRAWEVTHEADER");
0091   if (!m_mvtx_raw_event_header)
0092   {
0093     std::cout << PHWHERE << "::" << __func__ << ": Could not get MVTXRAWEVTHEADER from Node Tree" << std::endl;
0094     exit(1);
0095   }
0096   if (Verbosity() > 2)
0097   {
0098     m_mvtx_raw_event_header->identify();
0099   }
0100
0101   m_mvtx_raw_hit_container = findNode::getClass<MvtxRawHitContainer>(topNode, "MVTXRAWHIT");
0102   if (!m_mvtx_raw_hit_container)
0103   {
0104     std::cout << PHWHERE << "::" << __func__ << ": Could not get MVTXRAWHIT from Node Tree" << std::endl;
0105     exit(1);
0106   }
0107   if (Verbosity() > 2)
0108   {
0109     m_mvtx_raw_hit_container->identify();
0110   }
0111
0112   return Fun4AllReturnCodes::EVENT_OK;
0113 }
0114
0115 int MvtxFakeHitRate::process_event(PHCompositeNode* topNode)
0116 {
0117   // get the nodes
0118   if (m_num_strobes % 100000 == 0 && Verbosity() > 0)
0119   {
0120     std::cout << "MvtxFakeHitRate::process_event - processing strobe number: " << m_num_strobes << std::endl;
0121   }
0122   get_nodes(topNode);
0123
0124   // get lls from the MVTX raw event header
0125   for (unsigned int ihit = 0; ihit < m_mvtx_raw_hit_container->get_nhits(); ihit++)
0126   {
0127     // get this hit
0128     auto *mvtx_hit = m_mvtx_raw_hit_container->get_hit(ihit);
0129     if (!mvtx_hit)
0130     {
0131       std::cout << PHWHERE << "::" << __func__ << ": Could not get MVTX hit from container. Hit index: " << ihit << std::endl;
0132       continue;
0133     }
0134     if (m_hot_pixel_mask->is_masked(mvtx_hit))
0135     {
0136       m_masked_pixels_in_file = true;
0137     }
0138     // get the hit info
0139     uint64_t strobe = mvtx_hit->get_bco();
0140     uint8_t layer = mvtx_hit->get_layer_id();
0141     if (strobe > m_last_strobe)
0142     {
0143       m_last_strobe = strobe;
0144       m_num_strobes++;
0145     }
0146
0147     // if we are only looking at a specific layer, skip the rest
0148     if ((m_target_layer >= 0 && m_target_layer < 3) && (layer != m_target_layer))
0149     {
0150       continue;
0151     }
0152
0153     // generate pixel key
0154     MvtxPixelDefs::pixelkey this_pixelkey = MvtxPixelDefs::gen_pixelkey(mvtx_hit);
0155
0156     // add the hit to the hit map
0157     m_hit_map->add_hit(this_pixelkey, 1);
0158   }
0159
0160   return Fun4AllReturnCodes::EVENT_OK;
0161 }
0162
0163 int MvtxFakeHitRate::FillCurrentMaskTree()
0164 {
0165   // Fill the current mask tree
0166   m_current_masked_pixels.clear();
0167   m_current_nmasked = 0;
0168   std::vector<MvtxPixelDefs::pixelkey> hot_pixel_map = m_hot_pixel_mask->get_hot_pixel_map();
0169
0170   unsigned int nhits = m_hit_map->sum_hits(0);
0171   for (auto key : hot_pixel_map)
0172   {
0173     m_current_masked_pixels.push_back(key);
0174     m_current_nmasked++;
0175
0176     unsigned int nmasked_hits = m_hit_map->get_nhits(key);
0177     nhits -= nmasked_hits;
0178   }
0179
0180   m_current_threshold = calc_threshold(nhits);
0181   std::cout << "MvtxFakeHitRate::FillCurrentMaskTree - nmasked: " << m_current_nmasked << " threshold: " << m_current_threshold << std::endl;
0182
0183   m_current_mask->Fill();
0184
0185   return Fun4AllReturnCodes::EVENT_OK;
0186 }
0187
0188 double MvtxFakeHitRate::calc_threshold(int nhits) const
0189 {
0190   // Calculate the noise threshold
0191   if (nhits == 0)
0192   {
0193     return 0.0;
0194   }
0195   double npixels = 226492416.0;
0196   double denom = npixels * m_num_strobes;
0197   if (denom == 0)
0198   {
0199     return 0.0;
0200   }
0201   return (static_cast<double>(nhits) / denom);
0202 }
0203
0204 int MvtxFakeHitRate::CalcFHR()
0205 {
0206   // Calculate the fake hit rate
0207
0208   MvtxHitMap::pixel_hit_vector_t pixel_hit_vector = m_hit_map->get_pixel_hit_vector();
0209
0210   int npixels_with_hits = pixel_hit_vector.size();
0211   if (npixels_with_hits == 0)
0212   {
0213     std::cout << "MvtxFakeHitRate::CalcFHR - No pixels with hits" << std::endl;
0214     return Fun4AllReturnCodes::EVENT_OK;
0215   }
0216   m_max_masked_pixels = std::min(npixels_with_hits, m_max_masked_pixels);
0217
0218   // sort the pixel hit vector by hit count
0219   std::sort(pixel_hit_vector.begin(), pixel_hit_vector.end(), [](const MvtxHitMap::pixel_hits_pair_t& a, const MvtxHitMap::pixel_hits_pair_t& b)
0220             { return a.second > b.second; });
0221
0222   // get initial values
0223   unsigned int nhits_all = 0;
0224   m_masked_pixels.clear();
0225   for (auto & it : pixel_hit_vector)
0226   {
0227     nhits_all += it.second;
0228   }
0229   m_noise_threshold = calc_threshold(nhits_all);
0230   m_threshold_vs_nmasked->Fill(0.5, m_noise_threshold);
0231   m_tree->Fill();
0232
0233   for (int nmasked = 0; nmasked < m_max_masked_pixels; nmasked++)
0234   {
0235     m_masked_pixels.clear();
0236     m_num_masked_pixels = 0;
0237     int ipixel = 0;
0238     int nhits = nhits_all;
0239
0240     for (auto & it : pixel_hit_vector)
0241     {
0242       if (ipixel < nmasked)
0243       {
0244         m_masked_pixels.push_back(it.first);
0245         m_num_masked_pixels++;
0246         nhits -= it.second;
0247       }
0248       ipixel++;
0249     }
0250     if (nhits < 0)
0251     {
0252       std::cout << "MvtxFakeHitRate::CalcFHR - Error: nhits < 0" << std::endl;
0253       nhits = 0;
0254     }
0255     m_noise_threshold = calc_threshold(nhits);
0256
0257     m_threshold_vs_nmasked->Fill(m_num_masked_pixels + 0.5, m_noise_threshold);
0258     if (m_num_masked_pixels % 1000 == 0)
0259     {
0260       std::cout << "MvtxFakeHitRate::CalcFHR - nmasked: " << m_num_masked_pixels << " threshold: " << m_noise_threshold << std::endl;
0261     }
0262     m_tree->Fill();
0263   }
0264
0265   return Fun4AllReturnCodes::EVENT_OK;
0266 }
0267
0268 int MvtxFakeHitRate::End(PHCompositeNode* /*topNode*/)
0269 {
0270   std::cout << "MvtxFakeHitRate::End - Number of strobes: " << m_num_strobes << std::endl;
0271   // Fill the current mask tree
0272   FillCurrentMaskTree();
0273   if (m_masked_pixels_in_file)
0274   {
0275     std::cout << "MvtxFakeHitRate::End - Masked pixels found in file" << std::endl;
0276   }
0277   else
0278   {
0279     std::cout << "MvtxFakeHitRate::End - No masked pixels found in file" << std::endl;
0280   }
0281
0282   // Calculate the fake hit rate
0283   CalcFHR();
0284
0285   // Write the output
0286   std::cout << "MvtxFakeHitRate::End - Writing output to " << m_outputfile << std::endl;
0287   PHTFileServer::cd(m_outputfile);
0288   m_tree->Write();
0289   m_current_mask->Write();
0290   m_threshold_vs_nmasked->GetXaxis()->SetNdivisions(505);
0291   m_threshold_vs_nmasked->GetYaxis()->SetNdivisions(505);
0292   m_threshold_vs_nmasked->GetXaxis()->SetTitle("Masked pixels");
0293   m_threshold_vs_nmasked->GetYaxis()->SetTitle("(Active pix. / tot. pix.) per strobe");
0294   m_threshold_vs_nmasked->SetFillColor(kRed);
0295   m_threshold_vs_nmasked->SetLineColor(kRed);
0296   m_threshold_vs_nmasked->SetMarkerColor(kRed);
0297   m_threshold_vs_nmasked->Draw("PE1");
0298   m_threshold_vs_nmasked->Write();
0299   return Fun4AllReturnCodes::EVENT_OK;
0300 }