g4simulation/g4mvtx/PHG4MvtxDigitizer.cc

0001 // This is the new trackbase container version
0002
0003 #include "PHG4MvtxDigitizer.h"
0004
0005 // Move to new storage containers
0006 #include <trackbase/TrkrDefs.h>
0007 #include <trackbase/TrkrHit.h>  // for TrkrHit
0008 #include <trackbase/TrkrHitSet.h>
0009 #include <trackbase/TrkrHitSetContainer.h>
0010 #include <trackbase/TrkrHitTruthAssoc.h>
0011
0012 #include <g4detectors/PHG4CylinderGeom.h>
0013 #include <g4detectors/PHG4CylinderGeomContainer.h>
0014
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016 #include <fun4all/SubsysReco.h>  // for SubsysReco
0017
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHNode.h>  // for PHNode
0020 #include <phool/PHNodeIterator.h>
0021 #include <phool/PHRandomSeed.h>
0022 #include <phool/getClass.h>
0023 #include <phool/phool.h>  // for PHWHERE
0024
0025 #include <gsl/gsl_rng.h>  // for gsl_rng_alloc
0026
0027 #include <algorithm>
0028 #include <cstdlib>  // for exit
0029 #include <iostream>
0030 #include <set>
0031
0032 PHG4MvtxDigitizer::PHG4MvtxDigitizer(const std::string &name)
0033   : SubsysReco(name)
0034   , _energy_threshold(0.95e-6)
0035   , RandomGenerator(gsl_rng_alloc(gsl_rng_mt19937))
0036 {
0037   unsigned int seed = PHRandomSeed();  // fixed seed is handled in this funtcion
0038   std::cout << Name() << " random seed: " << seed << std::endl;
0039
0040   gsl_rng_set(RandomGenerator, seed);
0041
0042   if (Verbosity() > 0)
0043   {
0044     std::cout << "Creating PHG4MvtxDigitizer with name = " << name << std::endl;
0045   }
0046 }
0047
0048 PHG4MvtxDigitizer::~PHG4MvtxDigitizer()
0049 {
0050   gsl_rng_free(RandomGenerator);
0051 }
0052
0053 int PHG4MvtxDigitizer::InitRun(PHCompositeNode *topNode)
0054 {
0055   //-------------
0056   // Add Hit Node
0057   //-------------
0058   PHNodeIterator iter(topNode);
0059
0060   // Looking for the DST node
0061   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0062   if (!dstNode)
0063   {
0064     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0065     return Fun4AllReturnCodes::ABORTRUN;
0066   }
0067
0068   CalculateMvtxLadderCellADCScale(topNode);
0069
0070   //----------------
0071   // Report Settings
0072   //----------------
0073
0074   if (Verbosity() > 0)
0075   {
0076     std::cout << "====================== PHG4MvtxDigitizer::InitRun() =====================" << std::endl;
0077     for (auto &miter : _max_adc)
0078     {
0079       std::cout << " Max ADC in Layer #" << miter.first << " = " << miter.second << std::endl;
0080     }
0081     for (auto &miter : _energy_scale)
0082     {
0083       std::cout << " Energy per ADC in Layer #" << miter.first << " = " << 1.0e6 * miter.second << " keV" << std::endl;
0084     }
0085     std::cout << "===========================================================================" << std::endl;
0086   }
0087
0088   return Fun4AllReturnCodes::EVENT_OK;
0089 }
0090
0091 int PHG4MvtxDigitizer::process_event(PHCompositeNode *topNode)
0092 {
0093   // This code now only does the Mvtx
0094   DigitizeMvtxLadderCells(topNode);
0095
0096   return Fun4AllReturnCodes::EVENT_OK;
0097 }
0098
0099 void PHG4MvtxDigitizer::CalculateMvtxLadderCellADCScale(PHCompositeNode *topNode)
0100 {
0101   // defaults to 8-bit ADC, short-axis MIP placed at 1/4 dynamic range
0102
0103   PHG4CylinderGeomContainer *geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
0104
0105   if (!geom_container)
0106   {
0107     return;
0108   }
0109
0110   if (Verbosity())
0111   {
0112     std::cout << "Found CYLINDERGEOM_MVTX node" << std::endl;
0113   }
0114
0115   PHG4CylinderGeomContainer::ConstRange layerrange = geom_container->get_begin_end();
0116   for (PHG4CylinderGeomContainer::ConstIterator layeriter = layerrange.first;
0117        layeriter != layerrange.second;
0118        ++layeriter)
0119   {
0120     int layer = layeriter->second->get_layer();
0121     float thickness = (layeriter->second)->get_pixel_thickness();
0122     float pitch = (layeriter->second)->get_pixel_x();
0123     float length = (layeriter->second)->get_pixel_z();
0124
0125     float minpath = pitch;
0126     minpath = std::min(length, minpath);
0127     minpath = std::min(thickness, minpath);
0128     float mip_e = 0.003876 * minpath;
0129
0130     if (Verbosity())
0131     {
0132       std::cout << "mip_e = " << mip_e << std::endl;
0133     }
0134
0135     if (!_max_adc.contains(layer))
0136     {
0137       // cppcheck-suppress stlFindInsert
0138       _max_adc[layer] = 255;
0139       _energy_scale[layer] = mip_e / 64;
0140     }
0141   }
0142
0143   return;
0144 }
0145
0146 void PHG4MvtxDigitizer::DigitizeMvtxLadderCells(PHCompositeNode *topNode)
0147 {
0148   //----------
0149   // Get Nodes
0150   //----------
0151
0152   // new containers
0153   //=============
0154   // Get the TrkrHitSetContainer node
0155   TrkrHitSetContainer *trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0156   if (!trkrhitsetcontainer)
0157   {
0158     std::cout << "Could not locate TRKR_HITSET node, quit! " << std::endl;
0159     exit(1);
0160   }
0161
0162   // Get the TrkrHitTruthAssoc node
0163   auto *hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0164
0165   // Digitization
0166
0167   // We want all hitsets for the Mvtx
0168   TrkrHitSetContainer::ConstRange hitset_range = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::mvtxId);
0169   for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range.first;
0170        hitset_iter != hitset_range.second;
0171        ++hitset_iter)
0172   {
0173     // we have an itrator to one TrkrHitSet for the mvtx from the trkrHitSetContainer
0174     // get the hitset key so we can find the layer
0175     TrkrDefs::hitsetkey hitsetkey = hitset_iter->first;
0176     int layer = TrkrDefs::getLayer(hitsetkey);
0177     if (Verbosity() > 1)
0178     {
0179       std::cout << "PHG4MvtxDigitizer: found hitset with key: " << hitsetkey << " in layer " << layer << std::endl;
0180     }
0181
0182     // get all of the hits from this hitset
0183     TrkrHitSet *hitset = hitset_iter->second;
0184     TrkrHitSet::ConstRange hit_range = hitset->getHits();
0185     std::set<TrkrDefs::hitkey> hits_rm;
0186     for (TrkrHitSet::ConstIterator hit_iter = hit_range.first;
0187          hit_iter != hit_range.second;
0188          ++hit_iter)
0189     {
0190       TrkrHit *hit = hit_iter->second;
0191
0192       // Convert the signal value to an ADC value and write that to the hit
0193       // unsigned int adc = hit->getEnergy() / (TrkrDefs::MvtxEnergyScaleup *_energy_scale[layer]);
0194       if (Verbosity() > 0)
0195       {
0196         std::cout << "    PHG4MvtxDigitizer: found hit with key: " << hit_iter->first << " and signal " << hit->getEnergy() / TrkrDefs::MvtxEnergyScaleup << " in layer " << layer << std::endl;
0197       }
0198       // Remove the hits with energy under threshold
0199       bool rm_hit = false;
0200       if ((hit->getEnergy() / TrkrDefs::MvtxEnergyScaleup) < _energy_threshold)
0201       {
0202         if (Verbosity() > 0)
0203         {
0204           std::cout << "         remove hit, below energy threshold of " << _energy_threshold << std::endl;
0205         }
0206         rm_hit = true;
0207       }
0208       unsigned short adc = (unsigned short) (hit->getEnergy() / (TrkrDefs::MvtxEnergyScaleup * _energy_scale[layer]));
0209       adc = std::min(adc, _max_adc[layer]);
0210       hit->setAdc(adc);
0211
0212       if (rm_hit)
0213       {
0214         hits_rm.insert(hit_iter->first);
0215       }
0216     }
0217
0218     for (const auto &key : hits_rm)
0219     {
0220       if (Verbosity() > 0)
0221       {
0222         std::cout << "    PHG4MvtxDigitizer: remove hit with key: " << key << std::endl;
0223       }
0224       hitset->removeHit(key);
0225       if (hittruthassoc)
0226       {
0227         hittruthassoc->removeAssoc(hitsetkey, key);
0228       }
0229     }
0230   }
0231
0232   // end new containers
0233   //===============
0234
0235   return;
0236 }