Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
JEffectiveMassOffline1D.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3
4#include "TROOT.h"
5#include "TFile.h"
6#include "TH1D.h"
7
10
11#include "JAAnet/JVolume.hh"
14
15#include "JSupport/JSupport.hh"
18
19#include "JROOT/JManager.hh"
20
22
23#include "Jeep/JParser.hh"
24#include "Jeep/JMessage.hh"
25
26namespace {
27
28 const char* const Mass_t = "Mass";
29 const char* const Volume_t = "Volume";
30}
31
32/**
33 * \file
34 *
35 * Example program to histogram neutrino effective mass for triggered events for offline reconstruction files, using the event weights only.\n
36 * The unit of the contents of the histogram is \f$Mton\f$ or \f$km^{3}\f$.
37 * \author bjung, mdejong
38 */
39int main(int argc, char **argv)
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}
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
string outputFile
int main(int argc, char **argv)
Dynamic ROOT object management.
General purpose messaging.
#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
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Physics constants.
int numberOfBins
number of bins for average CDF integral of optical module
Definition JSirene.cc:73
ROOT TTree parameter settings of various packages.
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
Double_t getXmax() const
Get maximal abscissa value.
Definition JVolume.hh:102
Double_t getXmin() const
Get minimal abscissa value.
Definition JVolume.hh:91
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