diff --git a/src/programs/Simulation/mcsmear/ECALSmearer.cc b/src/programs/Simulation/mcsmear/ECALSmearer.cc index 7c469cce7..67ee978b4 100644 --- a/src/programs/Simulation/mcsmear/ECALSmearer.cc +++ b/src/programs/Simulation/mcsmear/ECALSmearer.cc @@ -6,7 +6,7 @@ //----------- // ecal_config_t (constructor) //----------- -ecal_config_t::ecal_config_t(const std::shared_ptr& event) { +ecal_config_t::ecal_config_t(const std::shared_ptr& event, const DECALGeometry *ecalGeom) { // Default Parameters @@ -19,23 +19,28 @@ ecal_config_t::ecal_config_t(const std::shared_ptr& event) { // Energy deposition in Geant ECAL_EN_GP0 = 1.71216e-2; - ECAL_EN_GP1 = 1.55070e-2; + ECAL_EN_GP1 = 1.e-2; ECAL_EN_GP2 = 0.0; // Time smearing factor ECAL_TSIGMA = 0.4; - - + // Single block energy threshold (applied after smearing) - ECAL_BLOCK_THRESHOLD = 15.0*k_MeV; + ECAL_ADC_THRESHOLD = 107; + // Baseline fluctuation + PED_SIGMA = 0; - // Get values from CCDB + ADC_EN_SCALE = 0.5447e-3; - cout << "Get ECAL/mc_energy parameters from CCDB..." << endl; + INT_OVER_PEAK = 5.18; + + // Get values from CCDB + cout << "Get ECAL/mc_energy parameters from CCDB ..." << endl; + map ecalparms; - + if(DEvent::GetCalib(event, "ECAL/mc_energy", ecalparms)) { jerr << "Problem loading ECAL/mc_energy from CCDB!" << endl; } else { @@ -44,27 +49,89 @@ ecal_config_t::ecal_config_t(const std::shared_ptr& event) { ECAL_EN_P0 = ecalparms["ECAL_EN_P0"]; ECAL_EN_P1 = ecalparms["ECAL_EN_P1"]; ECAL_EN_P2 = ecalparms["ECAL_EN_P2"]; - + ECAL_EN_GP0 = ecalparms["ECAL_EN_GP0"]; ECAL_EN_GP1 = ecalparms["ECAL_EN_GP1"]; ECAL_EN_GP2 = ecalparms["ECAL_EN_GP2"]; } - cout<<"get ECAL/mc_time parameters from calibDB"< ecaltime; + if(DEvent::GetCalib(event, "ECAL/mc_time", ecaltime)) { jerr << "Problem loading ECAL/mc_time from CCDB!" << endl; } else { ECAL_TSIGMA = ecaltime["ECAL_TSIGMA"]; } + cout << "Get ECAL/digi_scales parameters from CCDB ..." << endl; + map scale_factors; + + if (DEvent::GetCalib(event, "/ECAL/digi_scales", scale_factors)) + jout << "Error loading /ECAL/digi_scales !" << endl; + if (scale_factors.find("ADC_EN_SCALE") != scale_factors.end()) + ADC_EN_SCALE = scale_factors["ADC_EN_SCALE"]; + else + jerr << "Unable to get ADC_EN_SCALE from /ECAL/digi_scales !" << endl; + + cout << "Get ECAL/mc_parms from CCDB ..." << endl; + map mcparms; + + if(DEvent::GetCalib(event, "ECAL/mc_parms", mcparms)) { + jerr << "Problem loading ECAL/mc_parms from CCDB!" << endl; + } else { + ECAL_ADC_THRESHOLD = mcparms["THRESHOLD"]; + INT_OVER_PEAK = mcparms["INTEGRAL_PEAK"]; + PED_SIGMA = mcparms["PED_SIGMA"]; + } + + int max_chan = DECALGeometry::kECALMaxChannels; + + cout << "Get ECAL/gains from CCDB ..." << endl; + vector ecal_gains_ch; + + if (DEvent::GetCalib(event, "/ECAL/gains", ecal_gains_ch)){ + jout << "DECALHit_factory: Error loading /ECAL/gains !" << endl; + for (int ch = 0; ch < max_chan; ch ++) GAINS.push_back(1.); + } + else { + for (int ch = 0; ch < static_cast(ecal_gains_ch.size()); ch++) { + GAINS.push_back(ecal_gains_ch[ch]); + } + } + + cout << "Get ECAL/pedestals from CCDB ..." << endl; + vector ecal_pedestals_ch; + + if (DEvent::GetCalib(event, "/ECAL/pedestals", ecal_pedestals_ch)){ + jout << "DECALHit_factory: Error loading /ECAL/pedestals !" << endl; + for (int ch = 0; ch < max_chan; ch ++) PEDESTALS.push_back(100.); + } + else { + for (int ch = 0; ch < static_cast(ecal_pedestals_ch.size()); ch++) { + PEDESTALS.push_back(ecal_pedestals_ch[ch]); + } + } + + cout << "Get ECAL/bad_block from CCDB ..." << endl; + vector ecal_bad_blocks_ch; + + if (DEvent::GetCalib(event, "/ECAL/bad_block", ecal_bad_blocks_ch)){ + jout << "DECALHit_factory: Error loading /ECAL/bad_block !" << endl; + for (int ch = 0; ch < max_chan; ch ++) BAD_BLOCKS.push_back(0.); + } + else { + for (int ch = 0; ch < static_cast(ecal_bad_blocks_ch.size()); ch++) { + BAD_BLOCKS.push_back(ecal_bad_blocks_ch[ch]); + } + } + } //----------- -// SmearEvemt +// SmearEvent //----------- void ECALSmearer::SmearEvent(hddm_s::HDDM *record){ @@ -79,8 +146,15 @@ void ECALSmearer::SmearEvent(hddm_s::HDDM *record){ // A.S. new calibration of the ECAL double E = titer->getE(); double t = titer->getT(); + + int column=iter->getColumn(); + int row=iter->getRow(); - E *= ecal_config->ECAL_EN_SCALE; + int chan = ecalGeom->channel(row, column); + + double en_scale_cor = ecal_config->ECAL_EN_SCALE; + + E *= en_scale_cor; if(config->SMEAR_HITS) { @@ -91,24 +165,41 @@ void ECALSmearer::SmearEvent(hddm_s::HDDM *record){ // Subtract intrinsic Geant resolution double de_e_geant = pow(ecal_config->ECAL_EN_GP0/sqrt(E),2) + pow(ecal_config->ECAL_EN_GP1/E,2); - double sig_res = sqrt(de_e_expect - de_e_geant); + double sig_res = 0; + + if((de_e_expect - de_e_geant) > 0) + sig_res = sqrt(de_e_expect - de_e_geant); if(sig_res > 0) E *= (1. + gDRandom.SampleGaussian(sig_res)); - t += gDRandom.SampleGaussian(ecal_config->ECAL_TSIGMA); - + t += gDRandom.SampleGaussian(ecal_config->ECAL_TSIGMA); } + // A.S. Calculate energy threshold + int bad_block = ecal_config->BAD_BLOCKS.at(chan); + double adc_threshold = ecal_config->ECAL_ADC_THRESHOLD; + double pedestal = ecal_config->PEDESTALS.at(chan); + double pedestal_sigma = ecal_config->PED_SIGMA; + double gain = ecal_config->GAINS.at(chan); + double adc_en_scale = ecal_config->ADC_EN_SCALE; + double intOverPeak = ecal_config->INT_OVER_PEAK; + double pedestal_fluct = gDRandom.SampleGaussian(pedestal_sigma); + + double threshold = adc_en_scale*intOverPeak*en_scale_cor*gain*(adc_threshold - pedestal + pedestal_fluct); - // A.S. Don't apply energy threshold at the moment + double baseline_shift = adc_threshold - pedestal; + + if(bad_block > 0) continue; // Bad block - // if (E > ecal_config->ECAL_BLOCK_THRESHOLD) { - hddm_s::EcalHitList hits = iter->addEcalHits(); - hits().setE(E); - hits().setT(t); - // } + // Check the energy threshold. Do not produce hits if the baseline is above the threshold. + if ((E > threshold) && (baseline_shift > 0)) { // Check block threshold and baseline + hddm_s::EcalHitList hits = iter->addEcalHits(); + hits().setE(E); + hits().setT(t); + + } } if (config->DROP_TRUTH_HITS) @@ -119,5 +210,5 @@ void ECALSmearer::SmearEvent(hddm_s::HDDM *record){ if (ecals.size() > 0) ecals().deleteEcalTruthShowers(); } - + } diff --git a/src/programs/Simulation/mcsmear/ECALSmearer.h b/src/programs/Simulation/mcsmear/ECALSmearer.h index 72e843c89..8822086ca 100644 --- a/src/programs/Simulation/mcsmear/ECALSmearer.h +++ b/src/programs/Simulation/mcsmear/ECALSmearer.h @@ -5,10 +5,13 @@ #include "Smearer.h" +#include + + class ecal_config_t { public: - ecal_config_t(const std::shared_ptr& event); + ecal_config_t(const std::shared_ptr& event, const DECALGeometry *ecalGeom); double ECAL_EN_SCALE; @@ -24,10 +27,18 @@ class ecal_config_t // Time smearing factor double ECAL_TSIGMA; - // Single block energy threshold (applied after smearing) - double ECAL_BLOCK_THRESHOLD; - + + double ADC_EN_SCALE; + double INT_OVER_PEAK; + double ECAL_ADC_THRESHOLD; + + vector GAINS; + vector PEDESTALS; + vector BAD_BLOCKS; + + double PED_SIGMA = 0; + }; @@ -35,17 +46,18 @@ class ECALSmearer : public Smearer { public: ECALSmearer(const std::shared_ptr& event, mcsmear_config_t *in_config) : Smearer(event, in_config) { - ecal_config = new ecal_config_t(event); - } - ~ECALSmearer() { + event->GetSingle(ecalGeom); + ecal_config = new ecal_config_t(event,ecalGeom); + } + ~ECALSmearer() { delete ecal_config; } - + void SmearEvent(hddm_s::HDDM *record); private: ecal_config_t *ecal_config; - + const DECALGeometry *ecalGeom; };