#include #include #include #include "km3net-dataformat/offline/Head.hh" #include "km3net-dataformat/offline/Evt.hh" #include "JTools/JGrid.hh" #include "JTools/JPolint.hh" #include "JAAnet/JParticleTypes.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JAAnet/JDiffuseFluxHelper.hh" #include "JAAnet/JHondaFluxInterpolator.hh" #include "JSupport/JSupport.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JEvtWeightFileScannerSet.hh" #include "JOscProb/JOscParameters.hh" #include "JOscProb/JOscProbInterpolator.hh" #include "JOscProb/JOscProbHelper.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" #include "JROOT/JManager.hh" #include "TFile.h" #include "TH2D.h" /** * \file * Example program for plotting the atmospheric neutrino fluxes given by an azimuth-averaged Honda flux table\n * and comparing them with the atmospheric neutrino flux from the KM3NeT flux library. * * \author bjjung */ int main(int argc, char **argv) { using namespace std; using namespace JPP; typedef JGrid JGrid_t; JMultipleFileScanner_t inputFiles; string outputFile; string HondaTable; string oscProbTable; JOscParameters oscParameters; JGrid_t log10Egrid; JGrid_t coszGrid; int debug; try { JParser<> zap; zap['f'] = make_field(inputFiles); zap['T'] = make_field(HondaTable); zap['o'] = make_field(outputFile) = "honda.root"; zap['P'] = make_field(oscProbTable) = JPARSER::initialised(); zap['@'] = make_field(oscParameters) = JPARSER::initialised(); zap['X'] = make_field(log10Egrid) = make_grid(100, -1.0, 4.0); zap['C'] = make_field(coszGrid) = make_grid(200, -1.0, 1.0); zap['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } // Create atmospheric neutrino flux function JDiffuseFluxHelper flux1; JDiffuseFluxHelper flux2; const JFluxAtmospheric atmFlux; const JHondaFluxInterpolator2D<> interpolator(HondaTable.c_str()); if (!oscProbTable.empty()) { const JOscProbInterpolator<> oscProb(oscProbTable.c_str(), oscParameters); const JOscFlux oscFlux1(atmFlux, oscProb); const JOscFlux oscFlux2(interpolator, oscProb); flux1.configure(oscFlux1); flux2.configure(oscFlux2); } else { flux1.configure(atmFlux); flux2.configure(interpolator); } const vector nuTypes { TRACK_TYPE_ANTINUMU, TRACK_TYPE_ANTINUE, TRACK_TYPE_NUE, TRACK_TYPE_NUMU }; JManager h0(new TH2D("h0[%]", NULL, log10Egrid.getSize(), log10Egrid.getXmin(), log10Egrid.getXmax(), coszGrid .getSize(), coszGrid .getXmin(), coszGrid .getXmax())); JManager h1(new TH2D("h1[%]", NULL, log10Egrid.getSize(), log10Egrid.getXmin(), log10Egrid.getXmax(), coszGrid .getSize(), coszGrid .getXmin(), coszGrid .getXmax())); DEBUG(LEFT(15) << "E [GeV]" << RIGHT(10) << "cosz" << RIGHT(25) << "flux1" << RIGHT(52) << "flux2 [GeV^-1 * m^-2 * sr^-1 * s^-1]" << endl); for (vector::const_iterator type = nuTypes.cbegin(); type != nuTypes.cend(); ++type) { for (int i = 0; i < log10Egrid.getSize(); ++i) { const double X = log10Egrid.getX(i); for (int j = 0; j < coszGrid.getSize(); ++j) { const double cosz = coszGrid.getX(j); const double F0 = flux1(*type, X, cosz); const double F1 = flux2(*type, X, cosz); DEBUG(SCIENTIFIC(15,3) << pow(10.0,X) << FIXED (10,3) << cosz << SCIENTIFIC(25,3) << F0 << SCIENTIFIC(25,3) << F1 << endl); ASSERT(isfinite(F1)); h0[*type]->SetBinContent(i+1, j+1, F0); h1[*type]->SetBinContent(i+1, j+1, F1); } } } TFile out(outputFile.c_str(), "recreate"); out << h0 << h1; out.Write(); out.Close(); return 0; }