/* 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 "src/map/MapCppTrackerRadiationNoiseModel/MapCppTrackerRadiationNoiseModel.hh"
#include
#include "src/common_cpp/API/PyWrapMapBase.hh"
namespace MAUS {
PyMODINIT_FUNC init_MapCppTrackerRadiationNoiseModel(void) {
PyWrapMapBase::PyWrapMapBaseModInit
("MapCppTrackerRadiationNoiseModel", "", "", "", "");
}
MapCppTrackerRadiationNoiseModel::MapCppTrackerRadiationNoiseModel()
: MapBase("MapCppTrackerRadiationNoiseModel") {
// Do nothing
}
MapCppTrackerRadiationNoiseModel::~MapCppTrackerRadiationNoiseModel() {
delete _random;
// Do nothing
}
void MapCppTrackerRadiationNoiseModel::_birth(const std::string& argJsonConfigDocument) {
// Pull out the global settings
if (!Globals::HasInstance()) {
GlobalsManager::InitialiseGlobals(argJsonConfigDocument);
}
Json::Value* json = Globals::GetConfigurationCards();
// _clusters_on = (*json)["SciFiClusterReconOn"].asBool();
_random = new TRandom3();
// Build the geometery helper instance
MiceModule* module = Globals::GetReconstructionMiceModules();
std::vector modules =
module->findModulesByPropertyString("SensitiveDetector", "SciFi");
_geometry_helper = SciFiGeometryHelper(modules);
_geometry_helper.Build();
_cavityBrightness = (*json)["SciFiRadiationNoise_CavityBrightness"].asDouble();
// Make a list of the cavities in the experment
std::vector cavityPositions;
std::vector cavity_modules =
module->findModulesByPropertyExists("string", "CavityMode");
for (std::vector::const_iterator it = cavity_modules.begin();
it != cavity_modules.end(); ++it) {
cavityPositions.push_back((*it)->globalPosition().z());
}
if (cavityPositions.size() == 0) {
std::cerr << "No cavities found. No RF noise will be modelled.\n";
} else {
// Cycle through planes, find the nearest cavity and apply 1/r^2 reduction to
// the nominal rate.
for (std::vector::const_iterator it = modules.begin();
it != modules.end(); ++it) {
int tracker = (*it)->propertyInt("Tracker");
int station = (*it)->propertyInt("Station");
int plane = (*it)->propertyInt("Plane");
int plane_id = (station-1)*3+plane+1;
if (tracker == 0) plane_id *= -1;
double plane_z = (*it)->globalPosition().z();
double nearestCavity = 1.0e300;
for (std::vector::iterator cav_it = cavityPositions.begin();
cav_it != cavityPositions.end(); ++cav_it) {
double separation = fabs((*cav_it) - plane_z);
if (separation < nearestCavity) {
nearestCavity = separation;
}
}
// 1/r^2 losses assuming spherical emission
double planeRadius = (*it)->propertyDouble("ActiveRadius");
_planeNoiseRates[plane_id] = _cavityBrightness*planeRadius*planeRadius /
(4.0*nearestCavity*nearestCavity);
}
}
// Need a more precise determination - this will do for now!
// for (int plane_it = -15; plane_it <= 15; ++plane_it) {
// _planeNoiseRates[plane_it] = 100.0;
// }
_npe_mean = 5.0;
_npe_rms = 4.0;
_centralFibre = 106.5;
}
void MapCppTrackerRadiationNoiseModel::_death() {
// Do nothing
}
void MapCppTrackerRadiationNoiseModel::_process(Data* data) const {
Spill& spill = *(data->GetSpill());
/* return if not physics spill */
if (spill.GetDaqEventType() != "physics_event")
return;
if (spill.GetReconEvents()) {
for (unsigned int k = 0; k < spill.GetReconEvents()->size(); k++) {
SciFiEvent *event = spill.GetReconEvents()->at(k)->GetSciFiEvent();
if (!event) {
continue;
}
int spillNum = spill.GetSpillNumber();
int eventNum = spill.GetReconEvents()->at(k)->GetPartEventNumber();
for (std::map::const_iterator plane_it = _planeNoiseRates.begin();
plane_it != _planeNoiseRates.end(); ++plane_it) {
int plane_id = plane_it->first;
int tracker = (plane_id < 0?0:1);
int temp_plane_id = abs(plane_id);
int station = ((temp_plane_id - 1) / 3) + 1;
int plane = (temp_plane_id - 1) % 3;
int number_digits = _getNumberDigits(plane_id);
for (int digit_num = 0; digit_num < number_digits; ++digit_num) {
int channel = _getRandomChannel();
double npe = fabs(_random->Gaus(_npe_mean, _npe_rms));
double time = 0.0;
SciFiDigit* digit = new SciFiDigit(spillNum, eventNum,
tracker, station, plane, channel, npe, time);
event->add_digit(digit);
}
}
}
} else {
std::cout << "No recon events found\n";
}
}
int MapCppTrackerRadiationNoiseModel::_getNumberDigits(int plane) const {
std::map::const_iterator result = _planeNoiseRates.find(plane);
if (result == _planeNoiseRates.end()) {
throw Exceptions::Exception(Exceptions::nonRecoverable,
"Looking for a plane that doesn't exist!",
"MapCppTrackerRadiationNoiseModel::_getNumberDigits(int plane)");
}
double rate = _planeNoiseRates.find(plane)->second;
return _random->Poisson(rate);
}
int MapCppTrackerRadiationNoiseModel::_getRandomChannel() const {
return static_cast(_random->Uniform(0.0, _centralFibre*2.0));
}
} // ~namespace MAUS