#include #include #include #include #include "TROOT.h" #include "TFile.h" #include "TH1D.h" #include "TH2D.h" #include "TF1.h" #include "TFitResult.h" #include "JROOT/JMinimizer.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JTripod.hh" #include "JDetector/JTransmitter.hh" #include "JDetector/JHydrophone.hh" #include "JDetector/JModule.hh" #include "JDetector/JModuleRouter.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JTreeScanner.hh" #include "JROOT/JRootToolkit.hh" #include "JROOT/JManager.hh" #include "JAcoustics/JSoundVelocity.hh" #include "JAcoustics/JGeometry.hh" #include "JAcoustics/JEmitter.hh" #include "JAcoustics/JAcousticsToolkit.hh" #include "JAcoustics/JHit.hh" #include "JAcoustics/JEvent.hh" #include "JAcoustics/JEvt.hh" #include "JAcoustics/JEvtToolkit.hh" #include "JAcoustics/JSupport.hh" #include "JAcoustics/JTransmission_t.hh" #include "Jeep/JContainer.hh" #include "Jeep/JProperties.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * * Example program to plot acoustic fit results. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; typedef JContainer< vector > tripods_container; typedef JContainer< vector > transmitters_container; typedef JContainer< vector > hydrophones_container; typedef JContainer< set > disable_container; JMultipleFileScanner inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string detectorFile; string outputFile; JSoundVelocity V = getSoundVelocity; // default sound velocity tripods_container tripods; // tripods transmitters_container transmitters; // transmitters hydrophones_container hydrophones; // hydrophones int id; // emitter identifier disable_container disable; // disable tansmissions bool revert; string option; int numberOfEntries = 21; int debug; try { JParser<> zap("Example program to plot acoustic fit results."); JProperties properties; properties.insert(gmake_property(numberOfEntries)); zap['f'] = make_field(inputFile, "input file (output of JKatoomba[.sh])"); zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['a'] = make_field(detectorFile); zap['o'] = make_field(outputFile) = "canberra.root"; zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised(); zap['T'] = make_field(tripods, "tripod data"); zap['Y'] = make_field(transmitters, "transmitter data") = JPARSER::initialised(); zap['H'] = make_field(hydrophones, "hydrophone data") = JPARSER::initialised(); zap['M'] = make_field(getMechanics, "mechanics data") = JPARSER::initialised(); zap['E'] = make_field(id, "emitter identifier (-1 = all)") = -1; zap['!'] = make_field(disable, "disable transmission") = JPARSER::initialised(); zap['r'] = make_field(revert, "revert disabling of transmissions"); zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "LS"; zap['@'] = make_field(properties) = JPARSER::initialised(); zap['d'] = make_field(debug) = 2; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; } JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } JHashMap receivers; JHashMap emitters; for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) { receivers[i->getID()] = i->getLocation(); } for (tripods_container::const_iterator i = tripods.begin(); i != tripods.end(); ++i) { emitters[i->getID()] = JEmitter(i->getID(), i->getUTMPosition() - detector.getUTMPosition()); } for (transmitters_container::const_iterator i = transmitters.begin(); i != transmitters.end(); ++i) { try { emitters[i->getID()] = JEmitter(i->getID(), i->getPosition() + detector.getModule(i->getLocation()).getPosition()); } catch(const exception&) {} // if no module available, discard transmitter } const JGeometry geometry(detector, hydrophones); V.set(detector.getUTMZ()); JManager H2(new TH1D("[%].toa", NULL, 1000, -1.0e-2, +1.0e-2)); for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) { for (JHashMap::const_iterator i = emitters.begin(); i != emitters.end(); ++i) { if (id == i->first || id == -1) { H2[JTransmission_t(i->first, module->getID())]; } } } typedef JTreeScanner JTreeScanner_t; JTreeScanner_t in(inputFile); JTreeScanner_t::iterator p = in.begin(); while (inputFile.hasNext()) { STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); const JEvt* evt = inputFile.next(); const JModel model = getModel(*evt); DEBUG("model" << endl << model << endl); for ( ; p != in.end() && p-> begin()->getToA() < evt->UNIXTimeStart - 0.5; ++p) {} JTreeScanner_t::iterator q = p; for ( ; q != in.end() && q->rbegin()->getToA() <= evt->UNIXTimeStop + 0.5; ++q) {} DEBUG("events " << distance(p, q) << endl); for (JTreeScanner_t::iterator i = p; i != q; ++i) { if (id == i->getID() || id == -1) { if (emitters.has(i->getID())) { const JEmitter& emitter = emitters[i->getID()]; for (JEvent::const_iterator hit = i->begin(); hit != i->end(); ++hit) { if (receivers.has(hit->getID()) && geometry.hasLocation(receivers[hit->getID()])) { if ((disable.count(JTransmission_t(i->getID(), hit->getID())) == 0 && disable.count(JTransmission_t(-1, hit->getID())) == 0) == !revert) { const JLocation& location = receivers[hit->getID()]; if (geometry .has(location.getString()) && model.string.has(location.getString())) { const JGEOMETRY::JString& string = geometry [location.getString()]; const JMODEL ::JString& parameters = model.string[location.getString()]; const JPosition3D position = string.getPosition(parameters, location.getFloor()); const double D = emitter.getDistance(position); const double Vi = V.getInverseVelocity(D, emitter.getZ(), position.getZ()); const double toa = hit->getToE() + D * Vi; H2[JTransmission_t(i->getID(), hit->getID())]->Fill(hit->getToA() - toa); } } } } } } } p = q; } STATUS(endl); const floor_range range = getRangeOfFloors(detector); JManager HA(new TH2D("[%].mean", NULL, getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5, range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5)); JManager HB(new TH2D("[%].sigma", NULL, getNumberOfStrings(detector) + 0, -0.5, getNumberOfStrings(detector) - 0.5, range.getLength() + 1, range.getLowerLimit() - 0.5, range.getUpperLimit() + 0.5)); for (Int_t i = 1; i <= HA->GetXaxis()->GetNbins(); ++i) { HA->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first)); HB->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(geometry.at(i-1).first)); } const JModuleRouter router(detector); TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))"); const vector Q = { 0.33, 0.55, 0.66 }; for (JManager::iterator i = H2.begin(); i != H2.end(); ++i) { TH1* h1 = i->second; const JLocation location = router.getModule(i->first.rx).getLocation(); if (h1->GetEntries() >= numberOfEntries) { vector X(Q.size(), 0.0); h1->GetQuantiles(Q.size(), X.data(), Q.data()); f1.SetParameter(0, h1->GetMaximum()); f1.SetParameter(1, X[1]); f1.SetParameter(2, 0.5*(X[2] - X[0])); f1.SetParError(0, 0.1); f1.SetParError(1, 1.0e-6); f1.SetParError(2, 1.0e-6); TFitResultPtr result = h1->Fit(&f1, option.c_str(), "same", X[1] - 5.0 * (X[2] - X[0]), X[1] + 5.0 * (X[2] - X[0])); if (result.Get() == NULL) { ERROR("Invalid TFitResultPtr " << h1->GetName() << endl); } else { HA[i->first.tx]->Fill(geometry.getIndex(location.getString()), location.getFloor(), f1.GetParameter(1) * 1.0e3); HB[i->first.tx]->Fill(geometry.getIndex(location.getString()), location.getFloor(), f1.GetParameter(2) * 1.0e3); } } else { HA[i->first.tx]->Fill(geometry.getIndex(location.getString()), location.getFloor(), numeric_limits::lowest()); } } TFile out(outputFile.c_str(), "recreate"); out << H2 << HA << HB; out.Write(); out.Close(); }