/* 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 "src/common_cpp/Utils/JsonWrapper.hh"
#include "Interface/Squeak.hh"
#include "Utils/Exception.hh"
#include "Interface/dataCards.hh"
#include "Config/MiceModule.hh"
#include "src/common_cpp/API/PyWrapMapBase.hh"
#include "src/map/MapCppKLCellHits/MapCppKLCellHits.hh"
namespace MAUS {
PyMODINIT_FUNC init_MapCppKLCellHits(void) {
PyWrapMapBase::PyWrapMapBaseModInit
("MapCppKLCellHits", "", "", "", "");
}
MapCppKLCellHits::MapCppKLCellHits() : MapBase("MapCppKLCellHits") {
}
void MapCppKLCellHits::_birth(const std::string& argJsonConfigDocument) {
// Check if the JSON document can be parsed, else return error only
_stationKeys.push_back("kl");
// JsonCpp setup
Json::Value configJSON;
configJSON = JsonWrapper::StringToJson(argJsonConfigDocument);
// get the geometry
if (!configJSON.isMember("reconstruction_geometry_filename"))
throw(Exception(Exception::recoverable,
"Could not find geometry file",
"MapCppKLCellHits::birth"));
std::string filename;
filename = configJSON["reconstruction_geometry_filename"].asString();
// get the kl geometry modules
geo_module = new MiceModule(filename);
kl_modules = geo_module->findModulesByPropertyString("SensitiveDetector", "KL");
kl_mother_modules = geo_module->findModulesByPropertyString("Region", "KLregion");
}
void MapCppKLCellHits::_death() {}
void MapCppKLCellHits::_process(MAUS::Data* data) const {
// Get spill, break if there's no DAQ data
Spill *spill = data->GetSpill();
if (spill->GetReconEvents() == NULL)
return;
if (spill->GetDaqEventType() != "physics_event")
return;
ReconEventPArray *events = spill->GetReconEvents();
int nPartEvents = events->size();
for (int xPE = 0; xPE < nPartEvents; xPE++) {
KLDigitArray *kl_digits = (events->at(xPE))->GetKLEvent()
->GetKLEventDigitPtr()->GetKLDigitArrayPtr();
KLCellHitArray *kl_cellHits = (events->at(xPE))->GetKLEvent()
->GetKLEventCellHitPtr()->GetKLCellHitArrayPtr();
this->makeCellHits(kl_cellHits, kl_digits);
}
}
void MapCppKLCellHits::makeCellHits(KLCellHitArray* kl_CellHits, KLDigitArray* kl_digits) const {
int n_digits = kl_digits->size();
if (n_digits == 0)
return;
// Create a map of all hited PMTs.
std::map xDigitPos;
std::map::iterator it;
for (int xDigit = 0; xDigit < n_digits; xDigit++) {
std::string xKeyStr = (*kl_digits)[xDigit].GetKlKey();
KLChannelKey xKey(xKeyStr);
// Add this PMT to the map.
xDigitPos[xKeyStr] = xDigit;
}
while ( xDigitPos.size() > 1 ) {
// Get the first element of the map and check if we have a hit
// at the opposite side of the slab.
it = xDigitPos.begin();
KLChannelKey xKey(it->first);
// Get the opposite PMT coded as string.
std::string xOppositPmtKey_str = xKey.GetOppositeSidePMTStr();
if (xDigitPos.find(xOppositPmtKey_str) != xDigitPos.end()) {
// Create Cell hit.
KLCellHit xTheCellHit;
if (xKey.pmt() == 0) {
this->fillCellHit(xTheCellHit, (*kl_digits)[it->second],
(*kl_digits)[xDigitPos[xOppositPmtKey_str]]);
}
if (xKey.pmt() == 1) {
this->fillCellHit(xTheCellHit, (*kl_digits)[xDigitPos[xOppositPmtKey_str]],
(*kl_digits)[it->second]);
}
kl_CellHits->push_back(xTheCellHit);
// Erase both used hits from the map.
xDigitPos.erase(it);
xDigitPos.erase(xOppositPmtKey_str);
} else {
// Erese this hit from the map.
xDigitPos.erase(it);
}
}
}
void MapCppKLCellHits::fillCellHit(KLCellHit &cellHit, KLDigit &xDigit0, KLDigit &xDigit1) const {
std::string xKeyStr = xDigit0.GetKlKey();
KLChannelKey xKey(xKeyStr);
cellHit.SetCell(xKey.cell());
cellHit.SetDetector(xKey.detector());
cellHit.SetPartEventNumber(xDigit0.GetPartEventNumber());
cellHit.SetPhysEventNumber(xDigit0.GetPhysEventNumber());
// cell global position
// find the geo module corresponding to this hit
const MiceModule* hit_module = NULL;
Hep3Vector cellGlobalPos;
Hep3Vector cellErrorPos;
for ( unsigned int jj = 0; !hit_module && jj < kl_modules.size(); ++jj ) {
if ( kl_modules[jj]->propertyExists("Cell", "int") &&
kl_modules[jj]->propertyInt("Cell") == xKey.cell()) {
// got it
hit_module = kl_modules[jj];
} // end check on module
} // end loop over kl_modules
if (hit_module) {
cellGlobalPos = hit_module->globalPosition();
cellErrorPos = hit_module->dimensions()/sqrt(12);
} else {
cellGlobalPos.setX(-9999999.);
cellGlobalPos.setY(-9999999.);
cellGlobalPos.setZ(-9999999.);
cellErrorPos.setX(-9999999.);
cellErrorPos.setY(-9999999.);
cellErrorPos.setZ(-9999999.);
}
// get the local (relative to KL) cell positions from the geometry
const MiceModule* mother_module = NULL;
Hep3Vector cellLocalPos;
if (kl_mother_modules.size()) {
for ( unsigned int jj = 0; !mother_module && jj < kl_mother_modules.size(); ++jj ) {
mother_module= kl_mother_modules[jj];
}
cellLocalPos = hit_module->relativePosition(mother_module);
} else {
cellLocalPos.setX(-9999999.);
cellLocalPos.setY(-9999999.);
cellLocalPos.setZ(-9999999.);
}
cellHit.SetGlobalPosX(cellGlobalPos.x());
cellHit.SetGlobalPosY(cellGlobalPos.y());
cellHit.SetGlobalPosZ(cellGlobalPos.z());
cellHit.SetLocalPosX(cellLocalPos.x());
cellHit.SetLocalPosY(cellLocalPos.y());
cellHit.SetLocalPosZ(cellLocalPos.z());
cellHit.SetErrorX(cellErrorPos.x());
cellHit.SetErrorY(cellErrorPos.y());
cellHit.SetErrorZ(cellErrorPos.z());
// Charge of the digit can be unset because of the Zero suppresion of the fADCs.
int xChargeDigit0 = xDigit0.GetChargeMm();
int xChargeDigit1 = xDigit1.GetChargeMm();
cellHit.SetCharge((xChargeDigit0 + xChargeDigit1)/2);
if ((xChargeDigit0 + xChargeDigit1) == 0) {
cellHit.SetChargeProduct(0);
cellHit.SetFlag(false);
} else {
cellHit.SetChargeProduct(2 * xChargeDigit0 * xChargeDigit1 /
(xChargeDigit0 + xChargeDigit1));
cellHit.SetFlag(true);
}
}
}