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