Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
JEffectiveMassOffline1D.cc File Reference

Example program to histogram neutrino effective mass for triggered events for offline reconstruction files, using the event weights only. More...

#include <string>
#include <iostream>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "JAAnet/JVolume.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JEvtWeightFileScannerSet.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JROOT/JManager.hh"
#include "JPhysics/JConstants.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example program to histogram neutrino effective mass for triggered events for offline reconstruction files, using the event weights only.


The unit of the contents of the histogram is $Mton$ or $km^{3}$.

Author
bjung, mdejong

Definition in file JEffectiveMassOffline1D.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 39 of file JEffectiveMassOffline1D.cc.

40{
41 using namespace std;
42 using namespace JPP;
43 using namespace KM3NETDAQ;
44
45 JMultipleFileScanner_t inputFiles;
46 string outputFile;
47 bool logx;
48 int numberOfBins;
49 std::string option;
50 int debug;
51
52
53 try {
54
55 JParser<> zap("Example program to histogram neutrino effective mass for triggered events.");
56
57 zap['f'] = make_field(inputFiles);
58 zap['o'] = make_field(outputFile) = "Meff.root";
59 zap['X'] = make_field(logx, "Use logarithm of energy");
60 zap['N'] = make_field(numberOfBins, "Number of bins in the energy range of the MC simulation") = 10;
61 zap['O'] = make_field(option, "Result option") = Mass_t, Volume_t;
62 zap['d'] = make_field(debug) = 2;
63
64 zap(argc, argv);
65 }
66 catch(const exception &error) {
67 FATAL(error.what() << endl);
68 }
69
70 try {
71
72 double Xmin = numeric_limits<double>::max();
73 double Xmax = numeric_limits<double>::lowest();
74
75 JEvtWeightFileScannerSet<> scanners(inputFiles);
76
77 for (JEvtWeightFileScannerSet<>::const_iterator scanner = scanners.cbegin(); scanner != scanners.cend(); ++scanner) {
78
79 const JVolume volume(scanner->getHeader(), logx);
80
81 if (volume.getXmin() < Xmin) { Xmin = volume.getXmin(); }
82 if (volume.getXmax() > Xmax) { Xmax = volume.getXmax(); }
83 }
84
85 JManager<int, TH1D> hm(new TH1D("hm[%]", option.c_str(), numberOfBins, Xmin, Xmax));
86
87 for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
88
89 const JHead header = scanner->getHeader();
90
91 NOTICE("Scanning file type " << scanner->getName() << endl);
92 DEBUG (header << endl);
93
94 const JEvtWeight& weighter = scanner->getEvtWeighter();
95
96 if (!header.is_valid(&JHead::genvol)) {
97 FATAL("Mising normalisation in header." << endl);
98 }
99
100 double Wnorm = header.genvol.numberOfEvents;
101
102 if (header.is_valid(&JHead::cut_nu))
103 Wnorm *= 2*PI * (header.cut_nu.cosT.getUpperLimit() - header.cut_nu.cosT.getLowerLimit());
104 else
105 Wnorm *= 4*PI;
106
107 while (scanner->hasNext()) {
108
109 STATUS("event: " << setw(10) << scanner->getCounter() << '\r'); DEBUG(endl);
110
111 const Evt* event = scanner->next();
112 const Trk& primary = get_primary(*event);
113
114 double x = logx ? log10(primary.E) : primary.E;
115 double y = 0.0;
116 double w = (getVolume(weighter, *event) /
117 (logx ? log(10.0) * primary.E * hm->GetBinWidth(1) : hm->GetBinWidth(1)) /
118 Wnorm);
119
120 if (option == Mass_t) {
121 y = w * DENSITY_SEA_WATER * 1e-6; // Mton
122 } else if (option == Volume_t) {
123 y = w * 1e-9; // km^3
124 }
125
126 hm[primary.type]->Fill(x, y);
127 }
128 STATUS(endl);
129 }
130
131 TFile out(outputFile.c_str(), "recreate");
132
133 out << hm;
134
135 out.Close();
136 }
137 catch(const JException& error) {
138 FATAL(error << endl);
139 }
140}
string outputFile
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define NOTICE(A)
Definition JMessage.hh:64
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
int numberOfBins
number of bins for average CDF integral of optical module
Definition JSirene.cc:73
Monte Carlo run header.
Definition JHead.hh:1236
JAANET::cut_nu cut_nu
Definition JHead.hh:1596
const JHead & getHeader() const
Get header.
Definition JHead.hh:1270
JAANET::genvol genvol
Definition JHead.hh:1600
bool is_valid(T JHead::*pd) const
Check validity of given data member in JHead.
Definition JHead.hh:1319
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
Definition JManager.hh:47
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
double getVolume(const JType< JEvtWeightGSeaGen > &type, const Evt &evt)
Get volume of given event according given weighter.
const Trk & get_primary(const Evt &evt)
Auxiliary function to retrieve the primary track of an event.
static const double PI
Mathematical constants.
static const double DENSITY_SEA_WATER
Fixed environment values.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
Abstract base class for event weighing.
Definition JEvtWeight.hh:31
Auxiliary class for histogramming of effective volume.
Definition JVolume.hh:29
JRange_t cosT
Cosine zenith angle range
Definition JHead.hh:420
double numberOfEvents
Number of events.
Definition JHead.hh:721
Primary particle.
Definition JHead.hh:1174
int type
Particle type.
Definition JHead.hh:1204
Auxiliary class for organising Monte Carlo file scanners associated with event weighters.
std::vector< filescanner_type >::const_iterator const_iterator
std::vector< filescanner_type >::iterator iterator
Auxiliary base class for list of file names.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15