Jpp  18.0.0-rc.4
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JEffectiveMass1D.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 
11 
12 #include "JDAQ/JDAQEventIO.hh"
13 
16 
17 #include "JAAnet/JAAnetToolkit.hh"
18 #include "JAAnet/JHeadToolkit.hh"
19 #include "JAAnet/JVolume.hh"
20 
23 #include "JSupport/JSupport.hh"
24 
25 #include "JDetector/JDetector.hh"
27 
28 #include "Jeep/JParser.hh"
29 #include "Jeep/JMessage.hh"
30 
31 
32 /**
33  * \file
34  *
35  * Example program to histogram neutrino effective mass for triggered events inside the can volume.\n
36  * For genie/gSeaGen events the effective mass for triggered events inside the instrumented volume is also histogrammed.\n
37  * The unit of the contents of the histogram is Mton.
38  * \author mdejong, bstrandberg
39  */
40 int main(int argc, char **argv)
41 {
42  using namespace std;
43  using namespace JPP;
44  using namespace KM3NETDAQ;
45 
46  //------------------------------------------------------------------------------
47  // parse command line arguments, check that header can be found in the input file
48  //------------------------------------------------------------------------------
49 
50  JTriggeredFileScanner<> inputFile;
51  string outputFile;
52  double wall;
53  bool logx;
54  int numberOfBins;
55  string detectorFile;
56  int debug;
57 
58 
59  try {
60 
61  JParser<> zap("Example program to histogram neutrino effective mass for triggered events.");
62 
63  zap['f'] = make_field(inputFile);
64  zap['o'] = make_field(outputFile) = "Meff.root";
65  zap['R'] = make_field(wall, "Addition margin around the volume in which the considered events must reside") = 0.0;
66  zap['X'] = make_field(logx, "Use logarithm of energy");
67  zap['N'] = make_field(numberOfBins, "Number of bins in the energy range of the MC simulation") = 10;
68  zap['d'] = make_field(debug) = 1;
69  zap['a'] = make_field(detectorFile, "Detector file: if not provided, trigger fraction is not calculated") = "";
70 
71  zap(argc, argv);
72  }
73  catch(const exception &error) {
74  FATAL(error.what() << endl);
75  }
76 
77  //-------------------------------------------------------------------------
78  // deal with the header and init the detector
79  //-------------------------------------------------------------------------
80 
81  Head head;
82 
83  try {
84  head = getHeader(inputFile);
85  }
86  catch(const JException& error) {
87  FATAL(error);
88  }
89 
91 
92  if (detectorFile != "") {
93 
94  try {
95  load(detectorFile, detector);
96  }
97  catch(const JException& error) {
98  FATAL(error);
99  }
100  }
101 
102  //------------------------------------------------------------------------------
103  // create the can volume from the dimensions stored in the header and instrumented volume from detector file
104  //------------------------------------------------------------------------------
105 
106  const JHead buffer(head);
107 
108  const bool genie = is_gseagen(buffer);
109 
110  // the can volume comes from the generator header and the coordinate system is consistent by construction
111  const JVolume volume(head, logx);
112  JCylinder3D canvol = get<JCylinder3D>(buffer);
113  canvol.addMargin(wall);
114 
115  // as the detector is constructed from Jpp detector file, the coordinate system needs to be synced
116  // to the one used in the event generator
117  detector -= JAANET::getPosition( get<Vec>(buffer) );
118  JCylinder3D instvol( detector.begin(), detector.end() );
119  instvol.addMargin(wall);
120 
121  NOTICE("JVolume1D: Instrumented volume dimensions (zmin, zmax, r): " << instvol.getZmin() << " " << instvol.getZmax() << " " << instvol.getRadius() << endl );
122 
123  //------------------------------------------------------------------------------
124  // initialise output
125  //------------------------------------------------------------------------------
126 
127  TFile out(outputFile.c_str(), "recreate");
128 
129  TH1D hM("hM", "effective mass (can) [kg]" , numberOfBins, volume.getXmin(), volume.getXmax());
130  TH1D hm("hm", "effective mass (instrumented) [kg]", numberOfBins, volume.getXmin(), volume.getXmax());
131  hM.Sumw2();
132  hm.Sumw2();
133 
134  //------------------------------------------------------------------------------
135  // loop over triggered events in the input file
136  //------------------------------------------------------------------------------
137 
138  NOTICE("JVolume1D: Scanning triggered events." << endl);
139  while (inputFile.hasNext()) {
140 
141  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
142 
143  JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();
144 
145  const Evt* event = ps;
146 
147  if (has_neutrino(*event)) {
148 
149  const Trk& neutrino = get_neutrino(*event);
150  const JPosition3D& vertex = getPosition(neutrino);
151  const double E = neutrino.E;
152 
153  // genie/gSeaGen events filled with weight 1, histograms normalised after the second loop below
154  if ( genie ) {
155 
156  // events are filled to hM if they are inside the can
157  hM.Fill(volume.getX(E), canvol.is_inside(vertex) ? 1.0 : 0.0);
158 
159  // events are filled to hm if they are inside the instrumented volume
160  hm.Fill(volume.getX(E), instvol.is_inside(vertex) ? 1.0 : 0.0);
161 
162  }
163  // for other simulated events effective volume calculated from event weights, normalisation not required
164  else {
165  hM.Fill(volume.getX(E), volume.getW(hM.GetXaxis(), E));
166  }
167 
168  }
169  else {
170  WARNING("JEffectiveMass1D: cannot find neutrino in triggered event " << inputFile.getCounter() );
171  }
172 
173  } // end loop over triggered events
174  STATUS(endl);
175 
176  //------------------------------------------------------------------------------
177  // loop over generated events in the input file; FATAL if generator events have not been stored
178  //------------------------------------------------------------------------------
179 
180  if ( genie ) {
181 
182  // this histogram counts the generated events inside the can volume
183  TH1D* hNM = (TH1D*) hM.Clone("hNM");
184  hNM->Reset();
185 
186  // this histogram counts the generated events inside the instrumented volume
187  TH1D* hNm = (TH1D*) hm.Clone("hNm");
188  hNm->Reset();
189 
190  NOTICE("JVolume1D: Scanning generated events." << endl);
191 
192  JMultipleFileScanner<Evt> in(inputFile);
193 
194  while ( in.hasNext() ) {
195 
196  STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
197 
198  const Evt* event = in.next();
199 
200  if (has_neutrino(*event)) {
201 
202  const Trk& neutrino = get_neutrino(*event);
203  const JPosition3D& vertex = getPosition(neutrino);
204  const double E = neutrino.E;
205 
206  hNM->Fill(volume.getX(E), canvol.is_inside(vertex) ? 1.0 : 0.0);
207  hNm->Fill(volume.getX(E), instvol.is_inside(vertex) ? 1.0 : 0.0);
208 
209  }
210  else {
211  WARNING("JVolume1D: cannot find neutrino in generated event " << inputFile.getCounter() );
212  }
213 
214  }
215  STATUS(endl);
216 
217  if ( in.getCounter() == 0 ) {
218  FATAL("JVolume1D: generated events not stored in the input file, JTriggerEfficiency should be run without option -O");
219  }
220 
221  //------------------------------------------------------------------------------
222  // convert 'generated' and 'triggered' histograms to effective mass
223  //------------------------------------------------------------------------------
224 
225  hM.Divide(hNM);
226  hM.Scale(canvol.getVolume() * DENSITY_SEA_WATER * 1e-6); // Mton
227  hm.Divide(hNm);
228  hm.Scale(instvol.getVolume() * DENSITY_SEA_WATER * 1e-6); // Mton
229 
230  delete hNM;
231  delete hNm;
232  }
233 
234  out.Write();
235  out.Close();
236 
237 }
Utility class to parse command line options.
Definition: JParser.hh:1514
General exception.
Definition: JException.hh:23
Double_t getXmin() const
Get minimal abscissa value.
Definition: JVolume.hh:91
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
then usage $script< input file >[option[primary[working directory]]] nWhere option can be E
Definition: JMuonPostfit.sh:40
bool is_gseagen(const JHead &header)
Check for generator.
Definition: JHeadToolkit.hh:61
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
#define STATUS(A)
Definition: JMessage.hh:63
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Detector data structure.
Definition: JDetector.hh:89
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
Double_t getXmax() const
Get maximal abscissa value.
Definition: JVolume.hh:102
Double_t getX(const Double_t E, double constrain=false) const
Get abscissa value.
Definition: JVolume.hh:115
static const double DENSITY_SEA_WATER
Fixed environment values.
string outputFile
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
Data structure for detector geometry and calibration.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
Auxiliary class for histogramming of effective volume.
Definition: JVolume.hh:29
Cylinder object.
Definition: JCylinder3D.hh:39
Detector file.
Definition: JHead.hh:226
virtual bool hasNext() override
Check availability of next element.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
Double_t getW(TAxis *axis, const Double_t E) const
Get bin width corrected energy spectrum dependent weight.
Definition: JVolume.hh:156
JPosition3D getPosition(const Vec &pos)
Get position.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
#define NOTICE(A)
Definition: JMessage.hh:64
void addMargin(const double D)
Add (safety) margin.
Definition: JCylinder3D.hh:172
General purpose messaging.
Monte Carlo run header.
Definition: JHead.hh:1221
The Head class reflects the header of Monte-Carlo event files, which consists of keys (also referred ...
Definition: Head.hh:65
#define FATAL(A)
Definition: JMessage.hh:67
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
General purpose class for object reading from a list of file names.
Utility class to parse command line options.
virtual const pointer_type & next() override
Get next element.
counter_type getCounter() const
Get counter.
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
do set_variable DETECTOR_TXT $WORKDIR detector
General purpose class for multiple pointers.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY JAcoustics sh $DETECTOR_ID source JAcousticsToolkit sh CHECK_EXIT_CODE typeset A EMITTERS get_tripods $WORKDIR tripod txt EMITTERS get_transmitters $WORKDIR transmitter txt EMITTERS for EMITTER in
Definition: JCanberra.sh:46
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
then echo WARNING
Definition: JTuneHV.sh:91
int numberOfBins
number of bins for average CDF integral of optical module
Definition: JSirene.cc:69
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