g4simulation/g4waveformsim/CaloWaveformSim.cc

0001
0002 #include "CaloWaveformSim.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005
0006 #include <g4main/PHG4Hit.h>
0007 #include <g4main/PHG4HitContainer.h>
0008 #include <g4main/PHG4HitDefs.h>  // for hit_idbits
0009 #include <g4main/PHG4Particle.h>
0010 #include <g4main/PHG4TruthInfoContainer.h>
0011
0012 #include <cdbobjects/CDBTTree.h>  // for CDBTTree
0013
0014 #include <ffamodules/CDBInterface.h>
0015
0016 #include <ffaobjects/EventHeader.h>
0017
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHRandomSeed.h>
0020 #include <phool/getClass.h>
0021
0022 #include <calobase/TowerInfo.h>
0023 #include <calobase/TowerInfoContainer.h>
0024 #include <calobase/TowerInfoContainerSimv2.h>
0025
0026 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
0027 #include <g4detectors/PHG4CylinderCellGeom_Spacalv1.h>
0028 #include <g4detectors/PHG4CylinderGeomContainer.h>
0029 #include <g4detectors/PHG4CylinderGeom_Spacalv1.h>  // for PHG4CylinderGeom_Spaca...
0030 #include <g4detectors/PHG4CylinderGeom_Spacalv3.h>
0031
0032 #include <TF1.h>
0033 #include <TFile.h>
0034 #include <TProfile.h>
0035 #include <TSystem.h>
0036 #include <TTree.h>
0037 #include <algorithm>
0038 #include <cassert>
0039 #include <sstream>
0040 #include <string>
0041
0042 double CaloWaveformSim::template_function(double *x, double *par)
0043 {
0044   Double_t v1 = par[0] * h_template->Interpolate(x[0] - par[1]) + par[2];
0045   return v1;
0046 }
0047
0048 CaloWaveformSim::CaloWaveformSim(const std::string &name)
0049   : SubsysReco(name)
0050 {
0051 }
0052
0053 CaloWaveformSim::~CaloWaveformSim()
0054 {
0055   gsl_rng_free(m_RandomGenerator);
0056 }
0057
0058 int CaloWaveformSim::InitRun(PHCompositeNode *topNode)
0059 {
0060   // Initialize random generator
0061   m_RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
0062   unsigned int seed = PHRandomSeed();  // fixed seed handled in PHRandomSeed()
0063   gsl_rng_set(m_RandomGenerator, seed);
0064
0065   // Load waveform template
0066   const char *calibroot = getenv("CALIBRATIONROOT");
0067   if (!calibroot)
0068   {
0069     std::cerr << "CaloWaveformSim::InitRun missing CALIBRATIONROOT" << std::endl;
0070     exit(1);
0071   }
0072   std::string templatefilename = std::string(calibroot) + "/CaloWaveSim/" + m_templatefile;
0073   TFile *ft = TFile::Open(templatefilename.c_str());
0074   assert(ft && ft->IsOpen());
0075   h_template = static_cast<TProfile *>(ft->Get("hpwaveform"));
0076
0077   // Determine run number
0078   EventHeader *evtHeader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0079   m_runNumber = evtHeader ? evtHeader->get_RunNumber() : -1;
0080   if (Verbosity() > 0)
0081   {
0082     std::cout << "CaloWaveformSim::InitRun Run Number: " << m_runNumber << std::endl;
0083   }
0084
0085   // Detector-specific setup
0086   std::string url;
0087   if (m_dettype == CaloTowerDefs::CEMC)
0088   {
0089     m_detector = "CEMC";
0090     encode_tower = TowerInfoDefs::encode_emcal;
0091     decode_tower = TowerInfoDefs::decode_emcal;
0092     m_sampling_fraction = 2e-02;
0093     m_nchannels = 24576;
0094   }
0095   else if (m_dettype == CaloTowerDefs::HCALIN)
0096   {
0097     m_detector = "HCALIN";
0098     encode_tower = TowerInfoDefs::encode_hcal;
0099     decode_tower = TowerInfoDefs::decode_hcal;
0100     m_sampling_fraction = 0.162166;
0101     m_nchannels = 1536;
0102   }
0103   else  // HCALOUT
0104   {
0105     m_detector = "HCALOUT";
0106     encode_tower = TowerInfoDefs::encode_hcal;
0107     decode_tower = TowerInfoDefs::decode_hcal;
0108     m_sampling_fraction = 3.38021e-02;
0109     m_nchannels = 1536;
0110   }
0111
0112   // Gain settings
0113   // nobody understands this construct, please keep in mind that other
0114   // people have to read this and figure out what it does
0115   //  m_gain = m_highgain ? ((m_detector == "CEMC") ? 16 : 32) : 1;
0116   // use this instead:
0117   if (!m_highgain)
0118   {
0119     m_gain = 1;
0120   }
0121   else if (m_detector == "CEMC")
0122   {
0123     m_gain = 16;
0124   }
0125   else
0126   {
0127     m_gain = 32;
0128   }
0129
0130   // Data energy calibration
0131   if (!m_overrideCalibName)
0132   {
0133     m_calibName = m_detector + "_calib_ADC_to_ETower";
0134   }
0135   if (!m_overrideFieldName)
0136   {
0137     m_fieldname = m_detector + "_calib_ADC_to_ETower";
0138   }
0139   url = m_giveDirectURL ? m_directURL : CDBInterface::instance()->getUrl(m_calibName);
0140   if (!url.empty())
0141   {
0142     cdbttree = new CDBTTree(url);
0143   }
0144   else
0145   {
0146     std::cerr << "CaloWaveformSim::InitRun No data calibration for " << m_calibName << std::endl;
0147     exit(1);
0148   }
0149
0150   // MC energy calibration (optional)
0151   if (!m_overrideMCCalibName)
0152   {
0153     m_MC_calibName = m_detector + "_MC_RECALIB";
0154   }
0155   if (!m_overrideMCFieldName)
0156   {
0157     m_MC_fieldname = m_detector + "_calib_ADC_to_ETower";
0158   }
0159   url = m_giveDirectURL_MC ? m_directURL_MC : CDBInterface::instance()->getUrl(m_MC_calibName);
0160   if (!url.empty())
0161   {
0162     cdbttree_MC = new CDBTTree(url);
0163   }
0164   else if (Verbosity() > 0)
0165   {
0166     std::cout << "CaloWaveformSim::InitRun No MC calibration for " << m_MC_calibName << std::endl;
0167   }
0168
0169   // Time calibration (data)
0170   if (!m_overrideTimeCalibName)
0171   {
0172     m_calibName_time = m_detector + "_meanTime";
0173   }
0174   if (m_giveDirectURL_time)
0175   {
0176     url = m_directURL_time;
0177   }
0178   else
0179   {
0180     url = CDBInterface::instance()->getUrl(m_calibName_time);
0181     if (url.empty())
0182     {
0183       if (m_dotimecalib)
0184       {
0185         std::cerr << "CaloWaveformSim::InitRun No time calibration for " << m_calibName_time << std::endl;
0186         exit(1);
0187       }
0188     }
0189   }
0190   if (m_dotimecalib)
0191   {
0192     cdbttree_time = new CDBTTree(url);
0193   }
0194   if (Verbosity() > 0 && m_dotimecalib)
0195   {
0196     std::cout << "CaloWaveformSim::InitRun Time calibration from " << url << std::endl;
0197   }
0198
0199   // Time calibration (MC)
0200   if (!m_overrideMCTimeCalibName)
0201   {
0202     m_MC_calibName_time = m_detector + "_MC_meanTime";
0203   }
0204   if (m_giveDirectURL_MC_time)
0205   {
0206     url = m_directURL_MC_time;
0207     cdbttree_MC_time = new CDBTTree(url);
0208   }
0209   else
0210   {
0211     url = CDBInterface::instance()->getUrl(m_MC_calibName_time);
0212     if (!url.empty())
0213     {
0214       cdbttree_MC_time = new CDBTTree(url);
0215     }
0216     else if (m_dotimecalib)
0217     {
0218       std::cerr << "CaloWaveformSim::InitRun No MC time calibration for " << m_MC_calibName_time << std::endl;
0219       exit(1);
0220     }
0221   }
0222
0223   // Prepare waveform buffers
0224   m_waveforms.assign(m_nchannels, std::vector<float>(m_nsamples));
0225
0226   // Create node tree and finish
0227   CreateNodeTree(topNode);
0228   return Fun4AllReturnCodes::EVENT_OK;
0229 }
0230
0231 //____________________________________________________________________________..
0232 int CaloWaveformSim::process_event(PHCompositeNode *topNode)
0233 {
0234   if (Verbosity() > 0)
0235   {
0236     std::cout << "CaloWaveformSim::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0237   }
0238   // maybe we really need to get the geometry node in in every event(otherwise layergeom become invalid when we get to the second file in the list?):
0239   if (m_dettype == CaloTowerDefs::CEMC)
0240   {
0241     PHG4CylinderGeomContainer *layergeo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_CEMC");
0242     if (!layergeo)
0243     {
0244       std::cout << PHWHERE << " CYLINDERGEOM_CEMC Node missing, doing nothing." << std::endl;
0245       gSystem->Exit(1);
0246       exit(1);
0247     }
0248     const PHG4CylinderGeom *layergeom_raw = layergeo->GetFirstLayerGeom();
0249     assert(layergeom_raw);
0250
0251     layergeom = dynamic_cast<const PHG4CylinderGeom_Spacalv3 *>(layergeom_raw);
0252     assert(layergeom);
0253
0254     PHG4CylinderCellGeomContainer *seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_CEMC");
0255     if (!seggeo)
0256     {
0257       std::cout << PHWHERE << " CYLINDERCELLGEOM_CEMC Node missing, doing nothing." << std::endl;
0258       gSystem->Exit(1);
0259       exit(1);
0260     }
0261     PHG4CylinderCellGeom *geo_raw = seggeo->GetFirstLayerCellGeom();
0262     geo = dynamic_cast<PHG4CylinderCellGeom_Spacalv1 *>(geo_raw);
0263   }
0264
0265   // initialize the waveform
0266   for (auto &waveform : m_waveforms)
0267   {
0268     for (auto &sample : waveform)
0269     {
0270       sample = 0.;
0271     }
0272   }
0273   // waveform TH1
0274   TF1 *f_fit = new TF1(
0275       "f_fit", [this](double *x, double *par)
0276       { return this->template_function(x, par); },
0277       0, m_nsamples, 3);
0278   f_fit->SetParameter(0, 1.0);
0279   f_fit->SetParameter(1, 0.0);
0280   float template_peak = f_fit->GetMaximumX();
0281   float shift_of_shift = m_timeshiftwidth * gsl_rng_uniform(m_RandomGenerator);
0282
0283   float _shiftval = m_peakpos + shift_of_shift - template_peak;
0284
0285   f_fit->SetParameters(1, _shiftval, 0);
0286
0287   // get G4Hits
0288   std::string nodename = "G4HIT_" + m_detector;
0289   PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
0290   if (!hits)
0291   {
0292     std::cout << PHWHERE << " " << nodename << " Node missing, doing nothing." << std::endl;
0293     gSystem->Exit(1);
0294     exit(1);
0295   }
0296
0297   std::map<unsigned int,float> tbt_smear;
0298
0299
0300   // loop over hits
0301   for (PHG4HitContainer::ConstIterator hititer = hits->getHits().first; hititer != hits->getHits().second; hititer++)
0302   {
0303     PHG4Hit *hit = hititer->second;
0304     if (hit->get_t(1) - hit->get_t(0) > m_deltaT)
0305     {
0306       continue;
0307     }
0308
0309     // timing cut
0310
0311     // get eta phi bin
0312     unsigned short etabin = 0;
0313     unsigned short phibin = 0;
0314     float correction = 1.;
0315     maphitetaphi(hit, etabin, phibin, correction);
0316     unsigned int key = encode_tower(etabin, phibin);
0317     float calibconst = cdbttree->GetFloatValue(key, m_fieldname);
0318     float e_vis = hit->get_light_yield();
0319     e_vis *= correction;
0320     if (m_smear_const)
0321     {
0322       auto it = tbt_smear.find(key);
0323
0324       if(it != tbt_smear.end())
0325       {
0326         e_vis *= it->second;
0327       }
0328       else
0329       {
0330         tbt_smear[key] =  1.0+ gsl_ran_gaussian(m_RandomGenerator,factor_const);
0331         e_vis *= tbt_smear[key];
0332       }
0333     }
0334     float e_dep = e_vis / m_sampling_fraction;
0335     float ADC = (calibconst != 0) ? e_dep / calibconst : 0.;
0336     ADC *= m_gain;
0337
0338     if (cdbttree_MC)
0339     {
0340       float MC_calibconst = cdbttree_MC->GetFloatValue(key, m_MC_fieldname);
0341       ADC *= MC_calibconst;
0342     }
0343
0344     // if we have a tower by tower mean time, we shift the simulated waveform peak to that accordingly
0345     if (m_dotimecalib)
0346     {
0347       float meantime = cdbttree_time->GetFloatValue(key, m_fieldname_time);
0348       float MCmeantime = cdbttree_MC_time->GetFloatValue(key, m_MC_fieldname_time);
0349       assert(m_peakpos == 6);  // the MC mean time is derived when m_peakpos is set to 6
0350       _shiftval = m_peakpos + shift_of_shift - template_peak + meantime - MCmeantime;
0351     }
0352
0353     float t0 = hit->get_t(0) / m_sampletime;
0354     unsigned int tower_index = decode_tower(key);
0355     // here I will add the truth matching part
0356     //  for the cell reco, the truth matching info relys on edep not light yield, I will be consistent here :)
0357     TowerInfo *tower = m_CaloWaveformContainer->get_tower_at_channel(tower_index);
0358     TowerInfo::EdepMap &edepMap = tower->get_hitEdepMap();
0359     TowerInfo::ShowerEdepMap &showerMap = tower->get_showerEdepMap();
0360
0361     int showerID = hit->get_shower_id();
0362     float hitEdep = hit->get_edep();
0363
0364     edepMap[hit->get_hit_id()] += hitEdep;
0365     showerMap[showerID] += hitEdep;
0366
0367     f_fit->SetParameters(ADC, _shiftval + t0, 0.);
0368     for (int i = 0; i < m_nsamples; i++)
0369     {
0370       m_waveforms.at(tower_index).at(i) += f_fit->Eval(i);
0371     }
0372   }
0373
0374   // do noise here and add to waveform
0375
0376   if (m_noiseType == NoiseType::NOISE_TREE)
0377   {
0378     std::string ped_nodename = "PEDESTAL_" + m_detector;
0379     m_PedestalContainer = findNode::getClass<TowerInfoContainer>(topNode, ped_nodename);
0380
0381     if (!m_PedestalContainer)
0382     {
0383       std::cout << PHWHERE << " " << ped_nodename << " Node missing, doing nothing." << std::endl;
0384       gSystem->Exit(1);
0385       exit(1);
0386     }
0387   }
0388
0389   for (int i = 0; i < m_nchannels; i++)
0390   {
0391     std::vector<float> m_waveform_pedestal;
0392     m_waveform_pedestal.resize(m_nsamples);
0393     if (m_noiseType == NoiseType::NOISE_TREE)
0394     {
0395       TowerInfo *pedestal_tower = m_PedestalContainer->get_tower_at_channel(i);
0396       float pedestal_mean = 0;
0397       for (int j = 0; j < m_nsamples; j++)
0398       {
0399         m_waveform_pedestal.at(j) = (j < m_pedestalsamples) ? pedestal_tower->get_waveform_value(j) : pedestal_tower->get_waveform_value(m_pedestalsamples - 1);
0400         pedestal_mean += m_waveform_pedestal.at(j);
0401       }
0402       pedestal_mean /= m_nsamples;
0403       for (int j = 0; j < m_nsamples; j++)
0404       {
0405         m_waveform_pedestal.at(j) = (m_waveform_pedestal.at(j) - pedestal_mean) * m_pedestal_scale + pedestal_mean;
0406       }
0407     }
0408     for (int j = 0; j < m_nsamples; j++)
0409     {
0410       if (m_noiseType == NoiseType::NOISE_TREE)
0411       {
0412         // TowerInfo *pedestal_tower = m_PedestalContainer->get_tower_at_channel(i);
0413         // m_waveforms.at(i).at(j) += (j < m_pedestalsamples) ? pedestal_tower->get_waveform_value(j) : pedestal_tower->get_waveform_value(m_pedestalsamples - 1);
0414         m_waveforms.at(i).at(j) += m_waveform_pedestal.at(j);
0415       }
0416       if (m_noiseType == NoiseType::NOISE_GAUSSIAN)
0417       {
0418         m_waveforms.at(i).at(j) += gsl_ran_gaussian(m_RandomGenerator, m_gaussian_noise);
0419       }
0420       if (m_noiseType == NoiseType::NOISE_NONE)
0421       {
0422         m_waveforms.at(i).at(j) += m_fixpedestal;
0423       }
0424       // saturate at 2^14 - 1
0425       m_waveforms.at(i).at(j) = std::min<__gnu_cxx::__alloc_traits<class std::allocator<float> >::value_type>(m_waveforms.at(i).at(j), 16383);
0426       m_waveforms.at(i).at(j) = std::max<__gnu_cxx::__alloc_traits<class std::allocator<float> >::value_type>(m_waveforms.at(i).at(j), 0);
0427
0428       m_CaloWaveformContainer->get_tower_at_channel(i)->set_waveform_value(j, m_waveforms.at(i).at(j));
0429     }
0430   }
0431   delete f_fit;
0432   return Fun4AllReturnCodes::EVENT_OK;
0433 }
0434
0435 void CaloWaveformSim::maphitetaphi(PHG4Hit *g4hit, unsigned short &etabin, unsigned short &phibin, float &correction)
0436 {
0437   if (m_dettype == CaloTowerDefs::CEMC)
0438   {
0439     int scint_id = g4hit->get_scint_id();
0440     PHG4CylinderGeom_Spacalv3::scint_id_coder decoder(scint_id);
0441     std::pair<int, int> tower_z_phi_ID = layergeom->get_tower_z_phi_ID(decoder.tower_ID, decoder.sector_ID);
0442     const int &tower_ID_z = tower_z_phi_ID.first;
0443     const int &tower_ID_phi = tower_z_phi_ID.second;
0444     PHG4CylinderGeom_Spacalv3::tower_map_t::const_iterator it_tower =
0445         layergeom->get_sector_tower_map().find(decoder.tower_ID);
0446     assert(it_tower != layergeom->get_sector_tower_map().end());
0447
0448     const int etabin_cell = geo->get_etabin_block(tower_ID_z);  // block eta bin
0449     const int sub_tower_ID_x = it_tower->second.get_sub_tower_ID_x(decoder.fiber_ID);
0450     const int sub_tower_ID_y = it_tower->second.get_sub_tower_ID_y(decoder.fiber_ID);
0451     unsigned short etabinshort = etabin_cell * layergeom->get_n_subtower_eta() + sub_tower_ID_y;
0452     unsigned short phibin_cell = tower_ID_phi * layergeom->get_n_subtower_phi() + sub_tower_ID_x;
0453     etabin = etabinshort;
0454     phibin = phibin_cell;
0455
0456     // correction for emcal fiber
0457     if (light_collection_model.use_fiber_model())
0458     {
0459       const double z = 0.5 * (g4hit->get_local_z(0) + g4hit->get_local_z(1));
0460       assert(!std::isnan(z));
0461       correction *= light_collection_model.get_fiber_transmission(z);
0462     }
0463
0464     if (light_collection_model.use_fiber_model())
0465     {
0466       const double x = it_tower->second.get_position_fraction_x_in_sub_tower(decoder.fiber_ID);
0467       const double y = it_tower->second.get_position_fraction_y_in_sub_tower(decoder.fiber_ID);
0468       correction *= light_collection_model.get_light_guide_efficiency(x, y);
0469     }
0470   }
0471   else if (m_dettype == CaloTowerDefs::HCALIN)
0472   {
0473     // int layer = (g4hit->get_hit_id() >> PHG4HitDefs::hit_idbits);
0474     unsigned int iphi = (unsigned int) (g4hit->get_hit_id() >> PHG4HitDefs::hit_idbits) / 4;
0475     unsigned int ieta = g4hit->get_scint_id();
0476
0477     etabin = ieta;
0478     phibin = iphi;
0479   }
0480   else if (m_dettype == CaloTowerDefs::HCALOUT)
0481   {
0482     // int layer = (g4hit->get_hit_id() >> PHG4HitDefs::hit_idbits);
0483     unsigned int iphi = (unsigned int) (g4hit->get_hit_id() >> PHG4HitDefs::hit_idbits) / 5;
0484     unsigned int ieta = g4hit->get_scint_id();
0485
0486     etabin = ieta;
0487     phibin = iphi;
0488   }
0489   else
0490   {
0491     std::cout << PHWHERE << " Invalid detector type " << m_dettype << std::endl;
0492     gSystem->Exit(1);
0493     exit(1);
0494   }
0495 }
0496
0497 //____________________________________________________________________________..
0498 int CaloWaveformSim::End(PHCompositeNode * /*topNode*/)
0499 {
0500   std::cout << "CaloWaveformSim::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0501   return Fun4AllReturnCodes::EVENT_OK;
0502 }
0503
0504 void CaloWaveformSim::CreateNodeTree(PHCompositeNode *topNode)
0505 {
0506   PHNodeIterator topNodeItr(topNode);
0507   // DST node
0508   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(topNodeItr.findFirst("PHCompositeNode", "DST"));
0509   if (!dstNode)
0510   {
0511     std::cout << "PHComposite node created: DST" << std::endl;
0512     dstNode = new PHCompositeNode("DST");
0513     topNode->addNode(dstNode);
0514   }
0515   PHNodeIterator nodeItr(dstNode);
0516   PHCompositeNode *DetNode;
0517   // enum CaloTowerDefs::DetectorSystem and TowerInfoContainer::DETECTOR are different!!!!
0518   TowerInfoContainer::DETECTOR DetectorEnum = TowerInfoContainer::DETECTOR::DETECTOR_INVALID;
0519   std::string DetectorNodeName;
0520
0521   if (m_dettype == CaloTowerDefs::CEMC)
0522   {
0523     DetectorEnum = TowerInfoContainer::DETECTOR::EMCAL;
0524     DetectorNodeName = "CEMC";
0525   }
0526   else if (m_dettype == CaloTowerDefs::HCALIN)
0527   {
0528     DetectorEnum = TowerInfoContainer::DETECTOR::HCAL;
0529     DetectorNodeName = "HCALIN";
0530   }
0531   else if (m_dettype == CaloTowerDefs::HCALOUT)
0532   {
0533     DetectorEnum = TowerInfoContainer::DETECTOR::HCAL;
0534     DetectorNodeName = "HCALOUT";
0535   }
0536   else
0537   {
0538     std::cout << PHWHERE << " Invalid detector type " << m_dettype << std::endl;
0539     gSystem->Exit(1);
0540     exit(1);
0541   }
0542   DetNode = dynamic_cast<PHCompositeNode *>(nodeItr.findFirst("PHCompositeNode", DetectorNodeName));
0543   if (!DetNode)
0544   {
0545     DetNode = new PHCompositeNode(DetectorNodeName);
0546     dstNode->addNode(DetNode);
0547   }
0548   m_CaloWaveformContainer = new TowerInfoContainerSimv2(DetectorEnum);
0549
0550   PHIODataNode<PHObject> *newTowerNode = new PHIODataNode<PHObject>(m_CaloWaveformContainer, "WAVEFORM_" + m_detector, "PHObject");
0551   DetNode->addNode(newTowerNode);
0552 }