#include #include #include #include #include #include #include "TROOT.h" #include "TFile.h" #include "km3net-dataformat/offline/Head.hh" #include "km3net-dataformat/offline/MultiHead.hh" #include "km3net-dataformat/offline/Evt.hh" #include "km3net-dataformat/offline/io_online.hh" #include "km3net-dataformat/definitions/trkmembers.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JPMTRouter.hh" #include "JDetector/JModuleRouter.hh" #include "JDetector/JPMTParametersMap.hh" #include "JDetector/JPMTAnalogueSignalProcessor.hh" #include "JDynamics/JDynamics.hh" #include "JCompass/JEvt.hh" #include "JCompass/JSupport.hh" #include "JAcoustics/JEvt.hh" #include "JAcoustics/JSupport.hh" #include "JTrigger/JHit.hh" #include "JTrigger/JHitToolkit.hh" #include "JAAnet/JHead.hh" #include "JAAnet/JHeadToolkit.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JDAQ/JDAQEventIO.hh" #include "JDAQ/JDAQSummarysliceIO.hh" #include "JSupport/JParallelFileScanner.hh" #include "JSupport/JTreeScanner.hh" #include "JSupport/JFileRecorder.hh" #include "JSupport/JMonteCarloFileSupportkit.hh" #include "JSupport/JSupportToolkit.hh" #include "JSupport/JSupport.hh" #include "JSupport/JMeta.hh" #include "JReconstruction/JEvt.hh" #include "JReconstruction/JEvtToolkit.hh" #include "Jeep/JeepToolkit.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * Auxiliary program to convert fit results to Evt format. * * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace KM3NETDAQ; typedef JTYPELIST::typelist calibration_types; typedef JMultipleFileScanner JCalibration_t; JMultipleFileScanner<> inputFile; JFileRecorder::typelist> outputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string detectorFile; JCalibration_t calibrationFile; double Tmax_s; JPMTParametersMap pmtParameters; JQualitySorter qualitySorter; size_t numberOfFits; int debug; try { JParser<> zap("Auxiliary program to convert fit results to Evt format.\ \nThe option -L corresponds to the name of a shared library \ \nand function so to rearrange the order of fit results."); zap['f'] = make_field(inputFile); zap['o'] = make_field(outputFile) = "aanet.root"; zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['a'] = make_field(detectorFile) = ""; zap['+'] = make_field(calibrationFile) = JPARSER::initialised(); zap['T'] = make_field(Tmax_s) = 100.0; zap['P'] = make_field(pmtParameters, "PMT simulation data (or corresponding file name)") = JPARSER::initialised(); zap['L'] = make_field(qualitySorter) = JPARSER::initialised(); zap['N'] = make_field(numberOfFits) = 0; zap['d'] = make_field(debug) = 2; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } if (detectorFile == "" && !calibrationFile.empty()) { FATAL("Missing detector file." << endl); } typedef JParallelFileScanner< JTypeList > JParallelFileScanner_t; typedef JParallelFileScanner_t::multi_pointer_type multi_pointer_type; JDetector detector; if (detectorFile != "") { try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } } unique_ptr dynamics; try { dynamics.reset(new JDynamics(detector, Tmax_s)); dynamics->load(calibrationFile); } catch(const exception& error) { if (!calibrationFile.empty()) { FATAL(error.what()); } } const JPMTRouter pmt_router(dynamics ? dynamics->getDetector() : detector); const JModuleRouter mod_router(dynamics ? dynamics->getDetector() : detector); outputFile.open(); Head header; try { header = getHeader(inputFile); } catch(const exception& error) {} JAANET::JHead buffer(header); // DAQ livetime is not yet set for real data, for Monte Carlo, it is set in JTriggerEfficiency if (buffer.pull(&JAANET::JHead::DAQ) == buffer.end()) { buffer.DAQ.livetime_s = getLivetime(inputFile.begin(), inputFile.end()); buffer.push(&JAANET::JHead::DAQ); copy(buffer, header); } if (detectorFile != "") { buffer.calibration.buffer = (dynamics ? calibration::dynamical() : calibration::statical()); buffer.push(&JAANET::JHead::calibration); copy(buffer, header); } outputFile.put(header); outputFile.put(JMeta(argc, argv)); counter_type number_of_events = 0; for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) { STATUS("Processing: " << *i << endl); JParallelFileScanner_t in(*i); JTreeScanner mc(*i); in.setLimit(inputFile.getLimit()); Vec center(0,0,0); int mc_run_id = 0; try { const JAANET::JHead head(getHeader(*i)); center = get(head); mc_run_id = head.start_run.run_id; } catch(const exception& error) {} while (in.hasNext()) { STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl); multi_pointer_type ps = in.next(); JDAQEvent* tev = ps; JFIT::JEvt* evt = ps; if (dynamics) { dynamics->update(*tev); // update detector } JFIT::JEvt::iterator __end = evt->end(); if (numberOfFits > 0) { advance(__end = evt->begin(), min(numberOfFits, evt->size())); } if (qualitySorter.is_valid()) { partial_sort(evt->begin(), __end, evt->end(), qualitySorter); } Evt out; // output event if (mc.getEntries() != 0) { Evt* event = mc.getEntry(tev->getCounter()); if (event != NULL) { out = *event; // import Monte Carlo true data if (out.mc_run_id == 0) { out.mc_run_id = mc_run_id; } for (vector::iterator i = out.mc_trks.begin(); i != out.mc_trks.end(); ++i) { i->pos += center; // offset true positions } } } read(out, *tev); // import DAQ data if (!pmt_router->empty()) { // import calibration data for (vector::iterator i = out.mc_hits.begin(); i != out.mc_hits.end(); ++i) { if (pmt_router.hasPMT(i->pmt_id)) { const JPMTAddress address = pmt_router.getAddress(i->pmt_id); const JPMTIdentifier id = pmt_router.getIdentifier(address); const JPMT& pmt = pmt_router.getPMT(address); i->dom_id = id.getID(); i->channel_id = id.getPMTAddress(); i->pos = getPosition (pmt); i->dir = getDirection(pmt); } else { FATAL("Missing PMT" << i->pmt_id << endl); } } } if (!mod_router->empty()) { // import calibration data for (vector::iterator i = out.hits.begin(); i != out.hits.end(); ++i) { if (mod_router.hasModule(i->dom_id)) { const JPMTIdentifier id(i->dom_id, i->channel_id); const JPMTParameters& parameters = pmtParameters.getPMTParameters(id); const JPMTAnalogueSignalProcessor cpu(parameters); const JPMT& pmt = mod_router.getPMT(id); JTRIGGER::JHit::setSlewing(parameters.slewing); const JTRIGGER::JHit hit(getTime(i->tdc, pmt.getCalibration()), i->tot); i->pmt_id = pmt.getID(); i->pos = getPosition (pmt); i->dir = getDirection(pmt); i->t = hit.getT(); i->tot = hit.getToT(); i->a = cpu.getNPE(i->tot); } else { FATAL("Missing module " << i->dom_id << endl); } } } struct : public set { inline int get_index(const JUUID& element) const { const_iterator i = this->find(element); if (i != this->end()) return std::distance(this->begin(), i); else return -1; } } uuid; for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) { uuid.insert(fit->getUUID()); } for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) { // import reconstruction data Trk trk; trk.id = uuid.get_index(fit->getUUID()); trk.pos = Vec(fit->getX(), fit->getY(), fit->getZ()); trk.dir = Vec(fit->getDX(), fit->getDY(), fit->getDZ()); trk.t = fit->getT(); trk.E = fit->getE(); trk.lik = fit->getQ(); trk.rec_type = JPP_RECONSTRUCTION_TYPE; trk.status = (fit->getStatus() == COMPLETE_CHAIN ? TRK_ST_FINALSTATE : TRK_ST_UNDEFINED); if (fit->hasParentUUID()) { trk.mother_id = uuid.get_index(fit->getParentUUID()); } for (JHistory::const_iterator i = fit->getHistory().begin(); i != fit->getHistory().end(); ++i) { trk.rec_stages.push_back(i->type); } for (int i = 0; i != fit->getN(); ++i) { trk.fitinf.push_back(fit->getW(i)); } trk.error_matrix = fit->getV(); out.trks.push_back(trk); } out.id = ++number_of_events; outputFile.put(out); } } STATUS(endl); JMultipleFileScanner::typelist> io(inputFile); io >> outputFile; outputFile.close(); }