#include #include #include #include "km3net-dataformat/offline/Evt.hh" #include "km3net-dataformat/offline/Trk.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JAAnet/JLorentzBoost.hh" #include "JAAnet/JEvtWeightToolkit.hh" #include "JOscProb/JOscParameters.hh" #include "JOscProb/JOscProbInterpolator.hh" #include "JSirene/JEventShapeVariables.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JSupport/JSupport.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JEvtWeightFileScannerSet.hh" #include "JSupport/JMonteCarloFileSupportkit.hh" #include "Jeep/JeepToolkit.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" #include "TFile.h" #include "TH1D.h" /** * \file * * Auxiliary program to plot event shape variables. * \author bjjung */ int main(int argc, char **argv) { using namespace std; using namespace JPP; JMultipleFileScanner_t inputFiles; string outputFile; string detectorFile; JLimit_t numberOfEvents; bool useBoost; bool useWeights; JFluxMap fluxMaps; string oscProbTable; JOscParameters<> oscParameters; int debug; try { JParser<> zap("Auxiliary program to test event shape variables."); zap['f'] = make_field(inputFiles); zap['o'] = make_field(outputFile); zap['a'] = make_field(detectorFile) = ""; zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['B'] = make_field(useBoost); zap['W'] = make_field(useWeights); zap['@'] = make_field(fluxMaps) = JPARSER::initialised(); zap['P'] = make_field(oscProbTable) = JPARSER::initialised(); zap['#'] = make_field(oscParameters) = JPARSER::initialised(); zap['d'] = make_field(debug) = 3; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } if (useWeights && numberOfEvents != JLimit::max()) { FATAL("Cannot apply weighting to limited number of events."); } // Load detector file JDetector detector; if (!detectorFile.empty()) { try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } } JCylinder3D fiducialVolume; if (!detector.empty()) { fiducialVolume.set(detector.begin(), detector.end()); } // Define histograms TFile out(outputFile.c_str(), "recreate"); TH1D hT("hT", "thrust", 50, 0.5, 1.0); TH1D hA("hA", "aplanarity", 50, 0.0, 0.5); TH1D hS("hS", "sphericity", 100, 0.0, 1.0); TH1D hc("hc", "circularity", 100, 0.0, 1.0); TH1D hp("hp", "planar flow", 100, 0.0, 1.0); TH1D hC("hC", "C-variable", 100, 0.0, 1.0); TH1D hD("hD", "D-variable", 100, 0.0, 1.0); TH1D hH10("hH10", "H10", 100, 0.0, 1.0); TH1D hH20("hH20", "H20", 100, 0.0, 1.0); TH1D hH30("hH30", "H30", 100, 0.0, 1.0); TH1D hH40("hH40", "H40", 100, 0.0, 1.0); TH1D hq("hq", "Logarithmic eigenvalue ratio hit inertia tensor", 200, -20.0, 0.0); // Create file scanners JEvtWeightFileScannerSet<> scanners(inputFiles); if (!oscProbTable.empty()) { JOscProbInterpolator<> interpolator(oscProbTable.c_str(), oscParameters); fluxMaps.configure(make_oscProbFunction(interpolator)); } if (scanners.setFlux(fluxMaps) == 0) { WARNING("No flux function set." << endl); } STATUS(RIGHT(15) << "header ID" << RIGHT(15) << "event" << RIGHT(15) << "weight" << endl); // Loop over events for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) { scanner->setLimit(numberOfEvents); const JHead& header = scanner->getHeader(); if (detector.empty()) { // Use can volume if detector file has not been specified fiducialVolume = get(header); } while (scanner->hasNext()) { Evt* event = scanner->next(); const double weight = (useWeights ? scanner->getWeight(*event) : 1.0); STATUS(RIGHT (15) << distance(scanners.begin(), scanner) << RIGHT (15) << scanner->getCounter() << SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl); if (useBoost) { boostToCOM(*event); } const JVertex3D vertex = getVertex(*event); const JSphericityTensor S1(*event, 1.0); const JSphericityTensor S2(*event, 2.0); const JFoxWolframMoments H(*event, fiducialVolume, 4); const JHitInertiaTensor Q(*event, vertex); const double thrust = getThrust(*event); const double aplanarity = S2.getAplanarity(); const double sphericity = S2.getSphericity(); const double circularity = S2.getCircularity(); const double planarflow = S2.getPlanarFlow(); const double C = S1.getC(); const double D = S1.getD(); const double q = Q.getEigenvalueRatio(); hT.Fill(thrust, weight); hA.Fill(aplanarity, weight); hS.Fill(sphericity, weight); hc.Fill(circularity, weight); hp.Fill(planarflow, weight); hC.Fill(C, weight); hD.Fill(D, weight); hH10.Fill(H[1]/H[0], weight); hH20.Fill(H[2]/H[0], weight); hH30.Fill(H[3]/H[0], weight); hH40.Fill(H[4]/H[0], weight); hq.Fill(log10(q), weight); } } STATUS(endl); out.Write(); out.Close(); return 0; }