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