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