From 09107f141a72db5028046dc3eb5565316549422f Mon Sep 17 00:00:00 2001 From: Sungwon Kim Date: Sat, 1 Jun 2024 16:22:43 +0900 Subject: [PATCH] Fix units in SD action plugin, add wavelength information in output file --- plugins/FiberDRCaloSDAction.cpp | 25 ++++++++++++++------- plugins/Geant4Output2EDM4hep_DRC.cpp | 33 ++++++++++++++++++++++++++-- 2 files changed, 48 insertions(+), 10 deletions(-) diff --git a/plugins/FiberDRCaloSDAction.cpp b/plugins/FiberDRCaloSDAction.cpp index c2999bd4d..4f67fcd31 100644 --- a/plugins/FiberDRCaloSDAction.cpp +++ b/plugins/FiberDRCaloSDAction.cpp @@ -55,12 +55,13 @@ namespace dd4hep { G4float fWavlenStep; G4float fTimeStep; - G4double wavToE(G4double wav) { return h_Planck*c_light/wav; } + G4double wavToE(G4double wav) { return CLHEP::h_Planck*CLHEP::c_light/wav; } int findWavBin(G4double en) { int i = 0; for ( ; i < fWavBin+1; i++) { - if ( en < wavToE( (fWavlenStart - static_cast(i)*fWavlenStep)*nm ) ) + // if ( en < wavToE( (fWavlenStart - static_cast(i)*fWavlenStep)*nm ) ) + if ( en < wavToE( (fWavlenStart - static_cast(i)*fWavlenStep)*CLHEP::nm ) ) break; } @@ -77,13 +78,15 @@ namespace dd4hep { return i; } - bool applyYellowFilter(G4double energy) { - if (energy >= 2.69531) return true; + bool applyYellowFilter(G4double en) { + double energy = (en * 1000000.); // MeV -> eV unit change + + if (energy >= 2.69531) return true; // Photon w/ E > 2.69531 eV filtered TF1* uniform = new TF1("f1", "pol0", 0, 1); // TODO :: Pointer -> change to unique pointer uniform->SetParameter(0, 1); double randVal = uniform->GetRandom(0,1); // TODO :: Change random engine to DDSim or G4 random! - const int N = 23; + const int N = 24; double graph_X[N] = { 1.37760, 1.45864, 1.54980, @@ -106,7 +109,8 @@ namespace dd4hep { 2.47968, 2.53029, 2.58300, - 2.63796 }; + 2.63796, + 2.69531 }; double graph_Y[N] = { 0.903 * 0.903, 0.903 * 0.903, 0.903 * 0.903, @@ -129,13 +133,16 @@ namespace dd4hep { 0.345 * 0.345, 0.207 * 0.207, 0.083 * 0.083, - 0.018 * 0.018 }; + 0.018 * 0.018, + 0. }; TGraph* FilterEffGraph = new TGraph(N, graph_X, graph_Y); // TODO :: Pointer -> change to unique pointer double FilterEff = FilterEffGraph->Eval((double) energy); delete uniform; delete FilterEffGraph; - + + // filter efficiency == probability of photon passing filter + // == Probability of random value (from uniform distribution with range 0 ~ 1) larger than filter efficiency return ( randVal > FilterEff ); } @@ -190,6 +197,7 @@ namespace dd4hep { dd4hep::DDSegmentation::VolumeID ceren = static_cast(m_segmentation->decoder()->get(cID, "c")); bool IsCeren = static_cast(ceren); if ( !(IsCeren) ) { + // std::cout << "Photon energy : " << energy << std::endl; if ( m_userData.applyYellowFilter(energy) ) return true; } @@ -206,6 +214,7 @@ namespace dd4hep { hit->photonCount(); int wavBin = m_userData.findWavBin(energy); + std::cout << "wavBin : " << wavBin << std::endl; hit->CountWavlenSpectrum(wavBin); int timeBin = m_userData.findTimeBin(hitTime); hit->CountTimeStruct(timeBin); diff --git a/plugins/Geant4Output2EDM4hep_DRC.cpp b/plugins/Geant4Output2EDM4hep_DRC.cpp index 3cb4f67d5..f9ec446f7 100644 --- a/plugins/Geant4Output2EDM4hep_DRC.cpp +++ b/plugins/Geant4Output2EDM4hep_DRC.cpp @@ -49,12 +49,14 @@ namespace dd4hep { using drcalopair_t = std::pair< edm4hep::RawCalorimeterHitCollection, edm4hep::RawTimeSeriesCollection> ; // Required info for IDEA DRC sim hit using drcalomap_t = std::map< std::string, drcalopair_t >; // Required info for IDEA DRC sim hit + using drcaloWavmap_t = std::map< std::string, edm4hep::RawTimeSeriesCollection >; std::unique_ptr m_file { }; podio::Frame m_frame { }; edm4hep::MCParticleCollection m_particles { }; trackermap_t m_trackerHits; calorimetermap_t m_calorimeterHits; drcalomap_t m_drcaloHits; + drcaloWavmap_t m_drcaloWaves; stringmap_t m_runHeader; stringmap_t m_eventParametersInt; stringmap_t m_eventParametersFloat; @@ -282,11 +284,15 @@ void Geant4Output2EDM4hep_DRC::commit( OutputContext& /* ctxt */) { m_frame.put( std::move(calorimeterHits.first), colName + "RawHit"); m_frame.put( std::move(calorimeterHits.second), colName + "TimeStruct"); } + for (auto it = m_drcaloWaves.begin(); it != m_drcaloWaves.end(); ++it) { + m_frame.put( std::move(it->second), it->first + "WaveLen"); + } m_file->writeFrame(m_frame, m_section_name); m_particles.clear(); m_trackerHits.clear(); m_calorimeterHits.clear(); m_drcaloHits.clear(); + m_drcaloWaves.clear(); m_frame = {}; return; } @@ -325,6 +331,7 @@ void Geant4Output2EDM4hep_DRC::begin(const G4Event* event) { m_trackerHits.clear(); m_calorimeterHits.clear(); m_drcaloHits.clear(); + m_drcaloWaves.clear(); } /// Data conversion interface for MC particles to EDM4hep format @@ -584,8 +591,9 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext& /*ctxt*/, Geant4Sensitive* sd = coll->sensitive(); int hit_creation_mode = sd->hitCreationMode(); // Create the hit container even if there are no entries! - auto& hits = m_calorimeterHits[colName]; - auto& DRhits = m_drcaloHits[colName]; + auto& hits = m_calorimeterHits[colName]; + auto& DRhits = m_drcaloHits[colName]; + auto& DRwaves = m_drcaloWaves[colName]; // DEBUG // printout(INFO,"Geant4Output2EDM4hep_DRC","+++ Saving EDM4hep : Using Geant4DRCalorimeter::Hit"); @@ -624,6 +632,7 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext& /*ctxt*/, // For DRC raw calo hit & raw time series auto rawCaloHits = DRhits.first->create(); auto rawTimeStruct = DRhits.second->create(); + auto rawWaveStruct = DRwaves->create(); float samplingT = hit->GetSamplingTime(); float timeStart = hit->GetTimeStart(); @@ -635,7 +644,18 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext& /*ctxt*/, rawTimeStruct.setCharge( static_cast(hit->GetPhotonCount()) ); rawTimeStruct.setCellID( hit->cellID ); + // abuse time series for saving wavelength spectrum (for R&D purpose) + float samplingW = hit->GetSamplingWavlen(); + float wavMax = hit->GetWavlenMax(); + float wavMin = hit->GetWavlenMin(); + auto& wavmap = hit->GetWavlenSpectrum(); + rawWaveStruct.setInterval(samplingW); + rawWaveStruct.setTime(wavMin); + rawWaveStruct.setCharge( static_cast(hit->GetPhotonCount()) ); + rawWaveStruct.setCellID( hit->cellID ); + unsigned nbinTime = static_cast((timeEnd-timeStart)/samplingT); + unsigned nbinWav = static_cast((wavMax-wavMin)/samplingW); int peakTime = 0.; int peakVal = 0; @@ -656,6 +676,15 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext& /*ctxt*/, rawTimeStruct.addToAdcCounts(count); } + for (unsigned iwav = 1; iwav < nbinWav+1; iwav++) { + int count = 0; + + if ( wavmap.find(iwav)!=wavmap.end() ) + count = wavmap.at(iwav); + + rawWaveStruct.addToAdcCounts(count); + } + rawCaloHits.setCellID( hit->cellID ); rawCaloHits.setAmplitude( hit->GetPhotonCount() ); rawCaloHits.setTimeStamp( peakTime-1 + static_cast(timeStart/samplingT) );