Jpp  18.6.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
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

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 }
Utility class to parse command line options.
Definition: JParser.hh:1711
General exception.
Definition: JException.hh:24
bool is_valid(T JHead::*pd) const
Check validity of given data member in JHead.
Definition: JHead.hh:1319
data_type w[N+1][M+1]
Definition: JPolint.hh:867
Double_t getXmin() const
Get minimal abscissa value.
Definition: JVolume.hh:91
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
double getVolume(const JType< JEvtWeightGSeaGen > &type, const Evt &evt)
Get volume of given event according given weighter.
#define STATUS(A)
Definition: JMessage.hh:63
double numberOfEvents
Number of events.
Definition: JHead.hh:721
const Trk & get_primary(const Evt &evt)
Get primary.
static const double DENSITY_SEA_WATER
Fixed environment values.
string outputFile
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
Auxiliary class for histogramming of effective volume.
Definition: JVolume.hh:29
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys...
Definition: JManager.hh:43
Abstract base class for event weighing.
Definition: JEvtWeight.hh:28
Monte Carlo run header.
Definition: JHead.hh:1234
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
set_variable E_E log10(E_{fit}/E_{#mu})"
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
#define NOTICE(A)
Definition: JMessage.hh:64
JAANET::cut_nu cut_nu
Definition: JHead.hh:1595
static const double PI
Mathematical constants.
JRange_t cosT
Cosine zenith angle range.
Definition: JHead.hh:420
JAANET::genvol genvol
Definition: JHead.hh:1599
#define FATAL(A)
Definition: JMessage.hh:67
int type
MC: particle type in PDG encoding.
Definition: Trk.hh:24
Auxiliary base class for list of file names.
std::vector< filescanner_type >::iterator iterator
Primary particle.
Definition: JHead.hh:1174
const JHead & getHeader() const
Get header.
Definition: JHead.hh:1270
Auxiliary class for organising Monte Carlo file scanners associated with event weighters.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
std::vector< filescanner_type >::const_iterator const_iterator
int numberOfBins
number of bins for average CDF integral of optical module
Definition: JSirene.cc:70
int debug
debug level
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62