g4simulation/g4intt/PHG4InttDigitizer.cc

0001 // This is the new trackbase container version
0002
0003 #include "PHG4InttDigitizer.h"
0004
0005 #include <intt/InttMapping.h>
0006
0007 #include <g4detectors/PHG4CylinderGeom.h>
0008 #include <g4detectors/PHG4CylinderGeomContainer.h>
0009
0010 // Move to new storage containers
0011 #include <trackbase/InttDefs.h>
0012 #include <trackbase/TrkrDefs.h>
0013 #include <trackbase/TrkrHit.h>  // for TrkrHit
0014 #include <trackbase/TrkrHitSet.h>
0015 #include <trackbase/TrkrHitSetContainer.h>
0016 #include <trackbase/TrkrHitTruthAssoc.h>
0017
0018 #include <phparameter/PHParameterInterface.h>  // for PHParameterInterface
0019
0020 #include <fun4all/Fun4AllBase.h>  // for Fun4AllBase::VERB...
0021 #include <fun4all/Fun4AllReturnCodes.h>
0022 #include <fun4all/SubsysReco.h>  // for SubsysReco
0023
0024 #include <phool/PHCompositeNode.h>
0025 #include <phool/PHNode.h>  // for PHNode
0026 #include <phool/PHNodeIterator.h>
0027 #include <phool/PHRandomSeed.h>
0028 #include <phool/getClass.h>
0029 #include <phool/phool.h>  // for PHWHERE
0030
0031 #include <TSystem.h>
0032
0033 #include <gsl/gsl_randist.h>
0034 #include <gsl/gsl_rng.h>  // for gsl_rng_alloc
0035
0036 #include <algorithm>
0037 #include <cassert>
0038 #include <cmath>
0039 #include <cstdlib>  // for exit
0040 #include <iostream>
0041 #include <set>
0042 #include <utility>
0043
0044 PHG4InttDigitizer::PHG4InttDigitizer(const std::string &name)
0045   : SubsysReco(name)
0046   , PHParameterInterface(name)
0047 {
0048   InitializeParameters();
0049   unsigned int seed = PHRandomSeed();  // fixed seed is handled in this funtcion
0050   std::cout << Name() << " random seed: " << seed << std::endl;
0051   RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
0052   gsl_rng_set(RandomGenerator, seed);
0053 }
0054
0055 PHG4InttDigitizer::~PHG4InttDigitizer()
0056 {
0057   gsl_rng_free(RandomGenerator);
0058 }
0059
0060 int PHG4InttDigitizer::InitRun(PHCompositeNode *topNode)
0061 {
0062   std::cout << "PHG4InttDigitizer::InitRun: detector = " << detector << std::endl;
0063
0064   //-------------
0065   // Add Hit Node
0066   //-------------
0067
0068   PHNodeIterator iter(topNode);
0069
0070   // Looking for the DST node
0071   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0072   if (!dstNode)
0073   {
0074     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0075     return Fun4AllReturnCodes::ABORTRUN;
0076   }
0077
0078   CalculateLadderCellADCScale(topNode);
0079
0080   // Create the run and par nodes
0081   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0082   PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
0083
0084   std::string paramnodename = "G4CELLPARAM_" + detector;
0085   std::string geonodename = "G4CELLGEO_" + detector;
0086
0087   UpdateParametersWithMacro();
0088   // save this to the run wise tree to store on DST
0089   PHNodeIterator runIter(runNode);
0090   PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", detector));
0091   if (!RunDetNode)
0092   {
0093     RunDetNode = new PHCompositeNode(detector);
0094     runNode->addNode(RunDetNode);
0095   }
0096   SaveToNodeTree(RunDetNode, paramnodename);
0097   // save this to the parNode for use
0098   PHNodeIterator parIter(parNode);
0099   PHCompositeNode *ParDetNode = dynamic_cast<PHCompositeNode *>(parIter.findFirst("PHCompositeNode", detector));
0100   if (!ParDetNode)
0101   {
0102     ParDetNode = new PHCompositeNode(detector);
0103     parNode->addNode(ParDetNode);
0104   }
0105   PutOnParNode(ParDetNode, geonodename);
0106   mNoiseMean = get_double_param("NoiseMean");
0107   mNoiseSigma = get_double_param("NoiseSigma");
0108   mEnergyPerPair = get_double_param("EnergyPerPair");
0109
0110   /// If user hasn't called with custom calibration, load default
0111   if (!m_badmap.OfflineLoaded())
0112   {
0113     m_badmap.Load(); // Method loads with default tag
0114   }
0115   if (Verbosity())
0116   {
0117     std::cout << "InttBadChannelMap size: " << m_badmap.size() << std::endl;
0118   }
0119   if (VERBOSITY_MORE <= Verbosity())
0120   {
0121     m_badmap.Print();
0122   }
0123
0124   //----------------
0125   // Report Settings
0126   //----------------
0127
0128   if (Verbosity() > 0)
0129   {
0130     std::cout << "====================== PHG4InttDigitizer::InitRun() =====================" << std::endl;
0131     std::cout << " Masking " << m_badmap.size() << " channels" << std::endl;
0132     for (auto &iter1 : _max_adc)
0133     {
0134       std::cout << " Max ADC in Layer #" << iter1.first << " = " << iter1.second << std::endl;
0135     }
0136     for (auto &iter2 : _energy_scale)
0137     {
0138       std::cout << " Energy per ADC in Layer #" << iter2.first << " = " << 1.0e6 * iter2.second << " keV" << std::endl;
0139     }
0140     std::cout << "===========================================================================" << std::endl;
0141   }
0142
0143   return Fun4AllReturnCodes::EVENT_OK;
0144 }
0145
0146 int PHG4InttDigitizer::process_event(PHCompositeNode *topNode)
0147 {
0148   DigitizeLadderCells(topNode);
0149
0150   return Fun4AllReturnCodes::EVENT_OK;
0151 }
0152
0153 void PHG4InttDigitizer::CalculateLadderCellADCScale(PHCompositeNode *topNode)
0154 {
0155   // FPHX 3-bit ADC, thresholds are set in "set_fphx_adc_scale".
0156
0157   // PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_INTT");
0158   PHG4CylinderGeomContainer *geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
0159
0160   // if (!geom_container || !cells) return;
0161   if (!geom_container)
0162   {
0163     return;
0164   }
0165
0166   PHG4CylinderGeomContainer::ConstRange layerrange = geom_container->get_begin_end();
0167   for (PHG4CylinderGeomContainer::ConstIterator layeriter = layerrange.first;
0168        layeriter != layerrange.second;
0169        ++layeriter)
0170   {
0171     int layer = layeriter->second->get_layer();
0172     if (_max_fphx_adc.find(layer) == _max_fphx_adc.end())
0173     {
0174       std::cout << "Error: _max_fphx_adc is not available." << std::endl;
0175       gSystem->Exit(1);
0176     }
0177     float thickness = (layeriter->second)->get_thickness();  // cm
0178     float mip_e = 0.003876 * thickness;                      // GeV
0179     _energy_scale.insert(std::make_pair(layer, mip_e));
0180   }
0181
0182   return;
0183 }
0184
0185 void PHG4InttDigitizer::DigitizeLadderCells(PHCompositeNode *topNode)
0186 {
0187   //---------------------------
0188   // Get common Nodes
0189   //---------------------------
0190
0191   // Get the TrkrHitSetContainer node
0192   TrkrHitSetContainer *trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0193   if (!trkrhitsetcontainer)
0194   {
0195     std::cout << "Could not locate TRKR_HITSET node, quit! " << std::endl;
0196     exit(1);
0197   }
0198
0199   // Get the TrkrHitTruthAssoc node
0200   auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0201
0202   //-------------
0203   // Digitization
0204   //-------------
0205
0206   // We want all hitsets for the Intt
0207   TrkrHitSetContainer::ConstRange hitset_range = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::inttId);
0208   for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range.first;
0209        hitset_iter != hitset_range.second;
0210        ++hitset_iter)
0211   {
0212     // we have an itrator to one TrkrHitSet for the intt from the trkrHitSetContainer
0213     // get the hitset key so we can find the layer
0214     TrkrDefs::hitsetkey hitsetkey = hitset_iter->first;
0215     const int layer = TrkrDefs::getLayer(hitsetkey);
0216     const int ladder_phi = InttDefs::getLadderPhiId(hitsetkey);
0217     const int ladder_z = InttDefs::getLadderZId(hitsetkey);
0218
0219     if (Verbosity() > 1)
0220     {
0221       std::cout << "PHG4InttDigitizer: found hitset with key: " << hitsetkey << " in layer " << layer << std::endl;
0222     }
0223     // get all of the hits from this hitset
0224     TrkrHitSet *hitset = hitset_iter->second;
0225     TrkrHitSet::ConstRange hit_range = hitset->getHits();
0226     std::set<TrkrDefs::hitkey> dead_hits;  // hits on dead channel
0227     for (TrkrHitSet::ConstIterator hit_iter = hit_range.first;
0228          hit_iter != hit_range.second;
0229          ++hit_iter)
0230     {
0231       ++m_nCells;
0232
0233       TrkrHit *hit = hit_iter->second;
0234       TrkrDefs::hitkey hitkey = hit_iter->first;
0235       int strip_col = InttDefs::getCol(hitkey);  // strip z index
0236       int strip_row = InttDefs::getRow(hitkey);  // strip phi index
0237
0238       InttNameSpace::Offline_s ofl;
0239       ofl.layer = layer;
0240       ofl.ladder_phi = ladder_phi;
0241       ofl.ladder_z = ladder_z;
0242       ofl.strip_x = strip_row;
0243       ofl.strip_y = strip_col;
0244
0245       // Apply deadmap here if desired
0246       if (m_badmap.OfflineLoaded() && m_badmap.IsBad(ofl))
0247       {
0248         ++m_nDeadCells;
0249
0250         if (Verbosity() >= VERBOSITY_MORE)
0251         {
0252           std::cout << PHWHERE << " - dead strip at layer " << layer << ": ";
0253           hit->identify();
0254         }
0255
0256         dead_hits.insert(hit_iter->first);  // store hitkey of dead channels to be remove later
0257         continue;
0258       }
0259
0260       if (_energy_scale.count(layer) > 1)
0261       {
0262         assert(!"Error: _energy_scale has two or more keys.");
0263       }
0264       const float mip_e = _energy_scale[layer];
0265
0266       //      std::vector<double>& vadcrange = _max_fphx_adc[layer];
0267       // for(int i=0;i<8;i++) std::cout<<"Digitizer:: vadcrange= "<<vadcrange[i]<<std::endl;
0268
0269       // c++ upper_bound finds the bin location above the test value (or vadcrange.end() if there isn't one)
0270       //      auto irange = std::upper_bound(vadcrange.begin(), vadcrange.end(),
0271       //          hit->getEnergy() / (TrkrDefs::InttEnergyScaleup * (double) mip_e));
0272       //          hit->getEnergy() / 100.0);
0273       //      int adc = (irange-vadcrange.begin())-1;
0274       //      if (adc == -1) adc = 0;
0275
0276       double k = 85.7 / (TrkrDefs::InttEnergyScaleup * (double) mip_e);
0277       double E = hit->getEnergy() * k;  // keV
0278
0279       double gain = 100.0;
0280       double offset = 280.0;
0281       double para = 1.0;
0282       double e_vol = (E * pow(10, 3) * 1.6 * pow(10, -19) * pow(10, 15) * gain / 3.6) + offset;
0283       double v_dac = para * (e_vol - 210.0) / 4.0;
0284
0285       if (v_dac < 30)
0286       {
0287         v_dac = 15;
0288       }
0289       else if (v_dac < 60)
0290       {
0291         v_dac = 30;
0292       }
0293       else if (v_dac < 90)
0294       {
0295         v_dac = 60;
0296       }
0297       else if (v_dac < 120)
0298       {
0299         v_dac = 90;
0300       }
0301       else if (v_dac < 150)
0302       {
0303         v_dac = 120;
0304       }
0305       else if (v_dac < 180)
0306       {
0307         v_dac = 150;
0308       }
0309       else if (v_dac < 210)
0310       {
0311         v_dac = 180;
0312       }
0313       else
0314       {
0315         v_dac = 210;
0316       }
0317
0318       hit->setAdc(v_dac);
0319       /*
0320             std::cout<<"Digitizer:: getEnergy = "<<hit->getEnergy()<<std::endl;
0321             std::cout<<"Digitizer:: Energy = "<<E<<std::endl;
0322             std::cout<<"Digitizer:: k = "<<k<<std::endl;
0323             //std::cout<<"Digitizer:: capa = "<<capa<<std::endl;
0324             std::cout<<"Digitizer:: E2V = "<<e_vol<<std::endl;
0325             std::cout<<"Digitizer:: V2DAC = "<<v_dac<<std::endl;
0326             //std::cout<<"Digitizer:: mip_e = "<<TrkrDefs::InttEnergyScaleup * (double) mip_e<<std::endl;
0327             //std::cout<<"Digitizer:: getE/mip_e = "<<hit->getEnergy() /(TrkrDefs::InttEnergyScaleup * (double) mip_e)<<std::endl;
0328             //std::cout<<"Digitizer:: InttEnergyScaleup = "<<TrkrDefs::InttEnergyScaleup<<std::endl;
0329       */
0330       // if(adc==0) adc=15;
0331       // else adc=adc*30;
0332       // hit->setAdc(adc);
0333       // std::cout<<"Digitizer:: adc= "<<adc<<std::endl;
0334
0335       if (Verbosity() > 2)
0336       {
0337         std::cout << "PHG4InttDigitizer: found hit with layer " << layer << " ladder_z " << ladder_z << " ladder_phi " << ladder_phi
0338                   << " strip_col " << strip_col << " strip_row " << strip_row << " adc " << hit->getAdc() << std::endl;
0339       }
0340     }  // end loop over hits in this hitset
0341
0342     // remove hits on dead channel in TRKR_HITSET and TRKR_HITTRUTHASSOC
0343     for (const auto &key : dead_hits)
0344     {
0345       if (Verbosity() > 2)
0346       {
0347         std::cout << " PHG4InttDigitizer: remove hit with key: " << key << std::endl;
0348       }
0349       hitset->removeHit(key);
0350
0351       if (hittruthassoc)
0352       {
0353         hittruthassoc->removeAssoc(hitsetkey, key);
0354       }
0355     }
0356   }  // end loop over hitsets
0357
0358   return;
0359 }
0360
0361 //! end of process
0362 int PHG4InttDigitizer::End(PHCompositeNode * /*topNode*/)
0363 {
0364   if (Verbosity() >= VERBOSITY_SOME)
0365   {
0366     std::cout << "PHG4InttDigitizer::End - processed "
0367               << m_nCells << " cell with "
0368               << m_nDeadCells << " dead cells masked"
0369               << " (" << 100. * m_nDeadCells / m_nCells << "%)" << std::endl;
0370   }
0371
0372   return Fun4AllReturnCodes::EVENT_OK;
0373 }
0374
0375 void PHG4InttDigitizer::SetDefaultParameters()
0376 {
0377   set_default_double_param("NoiseMean", 457.2);
0378   set_default_double_param("NoiseSigma", 166.6);
0379   set_default_double_param("EnergyPerPair", 3.62e-9);  // GeV/e-h
0380   return;
0381 }
0382
0383 float PHG4InttDigitizer::added_noise()
0384 {
0385   //  float noise = gsl_ran_gaussian(RandomGenerator, mNoiseSigma) + mNoiseMean;
0386   //  noise = (noise < 0) ? 0 : noise;
0387
0388   // Note the noise is bi-polar, i.e. can make ths signal fluctuate up and down.
0389   // Much of the mNoiseSigma as extracted in https://github.com/sPHENIX-Collaboration/coresoftware/pull/580
0390   // is statistical fluctuation from the limited calibration data. They does not directly apply here.
0391   float noise = gsl_ran_gaussian(RandomGenerator, mNoiseMean);
0392
0393   return noise;
0394 }
0395
0396 void PHG4InttDigitizer::set_adc_scale(const int &layer, std::vector<double> userrange)
0397 {
0398   if (userrange.size() != nadcbins)
0399   {
0400     std::cout << "Error: vector in set_fphx_adc_scale(vector) must have eight elements." << std::endl;
0401     gSystem->Exit(1);
0402   }
0403   std::sort(userrange.begin(), userrange.end());
0404   _max_fphx_adc.insert(std::make_pair(layer, userrange));
0405 }