Skip to content

Commit

Permalink
Fix units in SD action plugin, add wavelength information in output file
Browse files Browse the repository at this point in the history
  • Loading branch information
@swkim95
swkim95 committed Jun 1, 2024
1 parent f5ed2f2 commit 09107f1
Show file tree
Hide file tree
Showing 2 changed files with 48 additions and 10 deletions.
25 changes: 17 additions & 8 deletions plugins/FiberDRCaloSDAction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<float>(i)*fWavlenStep)*nm ) )
// if ( en < wavToE( (fWavlenStart - static_cast<float>(i)*fWavlenStep)*nm ) )
if ( en < wavToE( (fWavlenStart - static_cast<float>(i)*fWavlenStep)*CLHEP::nm ) )
break;
}

Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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 );
}

Expand Down Expand Up @@ -190,6 +197,7 @@ namespace dd4hep {
dd4hep::DDSegmentation::VolumeID ceren = static_cast<dd4hep::DDSegmentation::VolumeID>(m_segmentation->decoder()->get(cID, "c"));
bool IsCeren = static_cast<bool>(ceren);
if ( !(IsCeren) ) {
// std::cout << "Photon energy : " << energy << std::endl;
if ( m_userData.applyYellowFilter(energy) ) return true;
}

Expand All @@ -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);
Expand Down
33 changes: 31 additions & 2 deletions plugins/Geant4Output2EDM4hep_DRC.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<writer_t> 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;
Expand Down Expand Up @@ -282,11 +284,15 @@ void Geant4Output2EDM4hep_DRC::commit( OutputContext<G4Event>& /* 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;
}
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -584,8 +591,9 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext<G4Event>& /*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");
Expand Down Expand Up @@ -624,6 +632,7 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext<G4Event>& /*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();
Expand All @@ -635,7 +644,18 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext<G4Event>& /*ctxt*/,
rawTimeStruct.setCharge( static_cast<float>(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<float>(hit->GetPhotonCount()) );
rawWaveStruct.setCellID( hit->cellID );

unsigned nbinTime = static_cast<unsigned>((timeEnd-timeStart)/samplingT);
unsigned nbinWav = static_cast<unsigned>((wavMax-wavMin)/samplingW);
int peakTime = 0.;
int peakVal = 0;

Expand All @@ -656,6 +676,15 @@ void Geant4Output2EDM4hep_DRC::saveCollection(OutputContext<G4Event>& /*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<int>(timeStart/samplingT) );
Expand Down

0 comments on commit 09107f1

Please sign in to comment.