/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus * * MAUS is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * MAUS is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with MAUS. If not, see . * */ #include #include #include #include "TH1D.h" #include "TFile.h" #include "TRandom3.h" #include "src/common_cpp/API/PyWrapMapBase.hh" #include "src/map/MapCppTrackerMCCalibratedNoise/MapCppTrackerMCCalibratedNoise.hh" namespace MAUS { PyMODINIT_FUNC init_MapCppTrackerMCCalibratedNoise(void) { PyWrapMapBase::PyWrapMapBaseModInit ("MapCppTrackerMCCalibratedNoise", "", "", "", ""); } /***** Initialization *****/ MapCppTrackerMCCalibratedNoise::MapCppTrackerMCCalibratedNoise() : MapBase("MapCppTrackerMCCalibratedNoise"), _calibration(217*3*5*2) { } /***** Destructor *****/ MapCppTrackerMCCalibratedNoise::~MapCppTrackerMCCalibratedNoise() { } /***** Birth *****/ void MapCppTrackerMCCalibratedNoise::_birth(const std::string& argJsonConfigDocument) { if (!Globals::HasInstance()) { GlobalsManager::InitialiseGlobals(argJsonConfigDocument); } static MiceModule* _mice_modules = Globals::GetMonteCarloMiceModules(); SF_modules = _mice_modules->findModulesByPropertyString("SensitiveDetector", "SciFi"); _configJSON = Globals::GetConfigurationCards(); _calibrationNoiseFilename = (*_configJSON)["SciFi_CalibratedNoiseFile"].asString(); _npe_scaling_factor = (*_configJSON)["SciFi_CalibratedNoiseNPEScaling"].asDouble(); _rate_scaling_factor = (*_configJSON)["SciFi_CalibratedNoiseRateScaling"].asDouble(); std::string map_file = (*_configJSON)["SciFiMappingFileName"].asString(); std::string calib_file = (*_configJSON)["SciFiCalibrationFileName"].asString(); std::string scifi_config_dir = (*_configJSON)["SciFiConfigDir"].asString(); this->_load_calibrations(scifi_config_dir, calib_file, map_file); TFile infile(_calibrationNoiseFilename.c_str(), "READ"); if (infile.IsZombie()) { throw Exceptions::Exception(Exceptions::nonRecoverable, "Could not open calibration noise file: "+_calibrationNoiseFilename, "MapCppTrackerMCCalibratedNoise::_birth()"); } std::string tracker = "US"; for (int station = 1; station < 6; ++station) { for (int plane = 0; plane < 3; ++plane) { std::stringstream plot_name; plot_name << "pe_" << tracker << "_" << station << "_" << plane; TH1D* the_plot = static_cast(infile.Get(plot_name.str().c_str())); the_plot->SetDirectory(0); _calibrationNoiseData.push_back(the_plot); } } tracker = "DS"; for (int station = 1; station < 6; ++station) { for (int plane = 0; plane < 3; ++plane) { std::stringstream plot_name; plot_name << "pe_" << tracker << "_" << station << "_" << plane; TH1D* the_plot = static_cast(infile.Get(plot_name.str().c_str())); the_plot->SetDirectory(0); _calibrationNoiseData.push_back(the_plot); } } infile.Close(); } /***** Death *****/ void MapCppTrackerMCCalibratedNoise::_death() { } /***** Process *****/ void MapCppTrackerMCCalibratedNoise::_process(Data* data) const { Spill* spill = data->GetSpill(); int spillNumber = spill->GetSpillNumber(); if (spill->GetReconEvents()) { for (unsigned int k = 0; k < spill->GetReconEvents()->size(); k++) { SciFiEvent *event = spill->GetReconEvents()->at(k)->GetSciFiEvent(); if (!event) { continue; } SciFiDigitPArray digits = event->digits(); for (int tracker = 0; tracker < 2; ++tracker) { for (int station = 1; station < 6; ++station) { for (int plane = 0; plane < 3; ++plane) { int plotID = _calculatePlotID(tracker, station, plane); int NChannels = _calibrationNoiseData[plotID]->GetSize(); for (int channelID = 1; channelID <= NChannels; ++channelID) { double rand = static_cast(std::rand()) / static_cast(RAND_MAX); if (rand < _calibrationNoiseData[plotID]->GetBinContent(channelID)*_rate_scaling_factor) { std::pair result = _calculateNPE(tracker, station, plane, channelID); SciFiDigit *a_digit = new SciFiDigit(spillNumber, k, tracker, station, plane, \ channelID, result.first, 0.0); a_digit->set_adc(result.second); digits.push_back(a_digit); } } } } } event->set_digits(digits); } } } std::pair MapCppTrackerMCCalibratedNoise::_calculateNPE(int tracker, \ int station, int plane, int channel) const { static TRandom3 randomizer = TRandom3(); int chan_id = calc_chan_id(tracker, station, plane, channel); std::pair chaninfo = _calibration[chan_id]; double pedestal = chaninfo.first; double gain = chaninfo.second; if ( gain == 0 ) return std::make_pair(0.0, 0.0); double smearing = randomizer.Gaus(0.0, gain/5.0); double npe = 2 + static_cast(-std::log(randomizer.Rndm())/_npe_scaling_factor); double adc = pedestal + gain*npe + smearing; double new_npe = (adc - pedestal)/gain; return std::make_pair(new_npe, adc); } int MapCppTrackerMCCalibratedNoise::_calculatePlotID(int tracker, int station, int plane) const {\ return (tracker == 0 ? 0 : 15) + (station-1)*3 + plane; } bool MapCppTrackerMCCalibratedNoise::_load_calibrations(std::string confDir, \ std::string calib_file, std::string cabling_file) { std::string calib_name = confDir+"/files/calibration/"+calib_file; std::string cabling_name = confDir+"/files/cabling/"+cabling_file; std::ifstream cabling_inf(cabling_name.c_str()); if (!cabling_inf) { throw(Exceptions::Exception(Exceptions::recoverable, "Could not load Tracker Mapping.", "MapCppTrackerMCCalibratedNoise::load_mapping")); } std::string line; ChanMapLookup chan_map(64*128); int line_count = 0; while ( getline(cabling_inf, line) ) { ++line_count; std::istringstream ist1(line.c_str()); ChanMap cmap; ist1 >> cmap.board >> cmap.bank >> cmap.chan_ro >> cmap.tracker >> cmap.station >> cmap.plane >> cmap.channel >> cmap.extWG >> cmap.inWG >> cmap.WGfib; cmap.bank = cmap.board*4 + cmap.bank; int UId = calc_uid(cmap.chan_ro, cmap.bank); // Check the UId does not exceed the bounds of the vector holding the channel map data if (static_cast(UId) > (chan_map.size() - 1)) { continue; } chan_map[UId] = cmap; } std::ifstream calib_inf(calib_name.c_str()); if (!calib_inf) { throw(Exceptions::Exception(Exceptions::recoverable, \ "Could not load Tracker Calibration.", \ "MapCppTrackerMCCalibratedNoise::load_calibration")); } std::string calib((std::istreambuf_iterator(calib_inf)), \ std::istreambuf_iterator()); Json::Reader reader; Json::Value calibration_data; if (!reader.parse(calib, calibration_data)) return false; size_t n_channels = calibration_data.size(); for ( Json::Value::ArrayIndex i = 0; i < n_channels; ++i ) { int bank = calibration_data[i]["bank"].asInt(); int channel_n = calibration_data[i]["channel"].asInt(); double adc_pedestal = calibration_data[i]["adc_pedestal"].asDouble(); double adc_gain = calibration_data[i]["adc_gain"].asDouble(); int uid = calc_uid(channel_n, bank); int tracker = chan_map[uid].tracker; int station = chan_map[uid].station; int plane = chan_map[uid].plane; int channel = chan_map[uid].channel; if (tracker < 0) continue; _calibration[calc_chan_id(tracker, station, plane, channel)] = \ std::make_pair(adc_pedestal, adc_gain); } return true; } } // ~namespace MAUS