Jpp  19.1.0-rc.1
the software that should make you happy
Functions
JDataPostfit.cc File Reference

Program to histogram fit results from reconstructed data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TMath.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/online/JDAQModuleIdentifier.hh"
#include "km3net-dataformat/definitions/reconstruction.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JVolume.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JDAQ/JDAQEventIO.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

Program to histogram fit results from reconstructed data.

Author
bofearraigh, rgracia

Definition in file JDataPostfit.cc.

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 43 of file JDataPostfit.cc.

44 {
45  using namespace std;
46  using namespace JPP;
47 
48  using namespace KM3NETDAQ;
49 
50  typedef JTriggeredFileScanner<JEvt> JTriggeredFileScanner_t;
51  typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
52 
53  JTriggeredFileScanner<JEvt> inputFile;
54  JLimit_t& numberOfEvents = inputFile.getLimit();
55  string outputFile;
56  string option;
57  size_t numberOfPrefits; // The number of reconstructed tracks to plot per event (1 by default)
58  JQualitySorter quality_sorter; // The criterion by which to pick the "best" reconstructed track for each event.
59  int application;
60  JRange<double> energy_range;
61  JRange<double> coszenith_range;
62  int debug;
63 
64  try {
65 
66  JParser<> zap("Program to histogram fit results from reconstructed data.");
67 
68  zap['f'] = make_field(inputFile);
69  zap['o'] = make_field(outputFile) = "postfit.root";
70  zap['n'] = make_field(numberOfEvents) = JLimit::max();
71  zap['N'] = make_field(numberOfPrefits) = 1;
72  zap['A'] = make_field(application) = JMUONENERGY, JMUONPREFIT, JMUONSIMPLEX, JMUONGANDALF, JMUONSTART; //makes histograms after this stage of the reconstruction
73  zap['L'] = make_field(quality_sorter) = JPARSER::initialised();
74  zap['E'] = make_field(energy_range) = JRange<double>();
75  zap['c'] = make_field(coszenith_range) = JRange<double>(-1.0,1.0);
76  zap['d'] = make_field(debug) = 2;
77 
78  zap(argc, argv);
79  }
80 
81  catch(const JException& error) {
82  FATAL(error);
83  }
84 
85  const double rad_to_deg = 180/M_PI;
86  const double R_m = 700.0;
87  const int MAX_OVERLAYS = 300;
88  const double MAX_TRIGGERED_HITS = 2500;
89  const double E_RANGE_MIN = -1.0;
90  const double E_RANGE_MAX = 7;
91  const double Z_RANGE_MIN = coszenith_range.getLowerLimit();
92  const double Z_RANGE_MAX = coszenith_range.getUpperLimit();
93 
94  TFile out(outputFile.c_str(), "recreate");
95  TH1D job("job", NULL, 100, 0.5, 100.5);
96  TH1D hz ("hz", "; cos(#theta_{zenith})" , 21 , Z_RANGE_MIN, Z_RANGE_MAX);
97  TH1D ho ("ho", "; # of overlays" , MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
98  TH2D hzo ("hzo", "; cos(#theta_{zenith}) ; # of overlays" , 21, Z_RANGE_MIN, Z_RANGE_MAX, MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
99  TH2D hxy ("hxy", "; x position ; y position" , 201, -R_m, R_m, 201, -R_m, R_m);
100  TH1D hn ("hn", "; # of triggered hits" , 600, -0.5, MAX_TRIGGERED_HITS-0.5);
101  TH1D hN ("hN", "; JEnergy # of hits" , 600, -0.5, MAX_TRIGGERED_HITS-0.5);
102  TH1D hq ("hq", "; quality parameter" , 200, -200.0, 800);
103  TH1D hb0 ("hb0", "; log(beta0)" , 60, -1, 3.5);
104  TH1D he ("he", "; log(Ereco [GeV])" , 75, E_RANGE_MIN, E_RANGE_MAX);
105  TH2D heo ("heo", "; log(Ereco [GeV]) ; # of overlays" , 75, E_RANGE_MIN, E_RANGE_MAX, MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
106  TH2D hen ("hen", "; log(Ereco [GeV]) ; # of triggered hits" , 75, E_RANGE_MIN, E_RANGE_MAX, 600 , -0.5, MAX_TRIGGERED_HITS-0.5);
107  TH2D heN ("heN", "; log(Ereco [GeV]) ; JEnergy # of hits" , 75, E_RANGE_MIN, E_RANGE_MAX, 600 , -0.5, MAX_TRIGGERED_HITS-0.5);
108  TH2D hzq ("hzq", "; cos(#theta_{zenith}); quality" , 21, Z_RANGE_MIN, Z_RANGE_MAX, 600, -200.0, 800.0);
109  TH2D hzn ("hzn", "; cos(#theta_{zenith}); # of triggered hits", 21, Z_RANGE_MIN, Z_RANGE_MAX, 600 , -0.5, MAX_TRIGGERED_HITS-0.5);
110  TH2D hzN ("hzN", "; cos(#theta_{zenith}); JEnergy # of hits" , 21, Z_RANGE_MIN, Z_RANGE_MAX, 600, -0.5, MAX_TRIGGERED_HITS-0.5);
111  TH2D hze ("hze", "; cos(#theta_{zenith}); log(Ereco [GeV])" , 21, Z_RANGE_MIN, Z_RANGE_MAX, 75, E_RANGE_MIN, E_RANGE_MAX);
112  TH2D hzb0("hzb0", "; cos(#theta_{zenith}); log(beta0)" , 21, Z_RANGE_MIN, Z_RANGE_MAX, 60, E_RANGE_MIN, 3.5);
113 
114 
115  while (inputFile.hasNext()) { // Loop over each reconstructed event
116 
117  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); // if debug > set value, print out this line
118 
119  multi_pointer_type ps = inputFile.next();
120 
121  JDAQEvent* tev = ps;
122  JEvt* evt = ps;
123 
124  job.Fill(1.0);
125 
126  if (!evt->empty()) {
127 
128  JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_event(application));
129 
130  if (evt->begin() == __end) {
131  continue;
132  }
133 
134  if (numberOfPrefits > 0) {
135 
136  JEvt::iterator __q = __end;
137 
138  advance(__end = evt->begin(), min(numberOfPrefits, (size_t) distance(evt->begin(), __q)));
139 
140  partial_sort(evt->begin(), __end, __q, quality_sorter);
141 
142  } else {
143 
144  sort(evt->begin(), __end, quality_sorter);
145  }
146 
147  if (numberOfPrefits > 0) {
148  advance(__end = evt->begin(), min(numberOfPrefits, evt->size()));
149  }
150 
151  for (JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) {
152 
153  JTrack3D track = getTrack(*fit);
154 
155  if (track.getDZ() > Z_RANGE_MIN && track.getDZ() < Z_RANGE_MAX && fit->getE() > energy_range.getLowerLimit() && fit->getE() < energy_range.getUpperLimit() ) {
156  const double overlays = tev->getOverlays();
157  const double Efit = fit->getE();
158 
159  ho .Fill(overlays);
160  hn .Fill((double) tev->size<JDAQTriggeredHit>());
161  hz .Fill(track.getDZ());
162  hzn.Fill(track.getDZ(), (double) tev->size<JDAQTriggeredHit>());
163  hzo.Fill(track.getDZ(), overlays) ;
164  hxy.Fill(track.getX(), track.getY());
165  hq .Fill(fit->getQ());
166  hzq.Fill(track.getDZ(), fit->getQ() );
167 
168  if (fit->hasW(JENERGY_NUMBER_OF_HITS)) {
169  hN .Fill( fit->getW(JENERGY_NUMBER_OF_HITS) );
170  heN .Fill(log10(Efit), fit->getW(JENERGY_NUMBER_OF_HITS) );
171  hzN .Fill(track.getDZ(), fit->getW(JENERGY_NUMBER_OF_HITS) );
172  }
173 
174  if (fit->hasW(JGANDALF_BETA0_RAD)) {
175  hb0 .Fill(log10( rad_to_deg* fit->getW(JGANDALF_BETA0_RAD) ));
176  hzb0 .Fill(track.getDZ(), log10( rad_to_deg* fit->getW(JGANDALF_BETA0_RAD)) );
177  }
178  he .Fill(log10(Efit));
179  hen .Fill(log10(Efit), (double) tev->size<JDAQTriggeredHit>());
180  hze .Fill(track.getDZ(), log10(Efit) );
181  heo .Fill(log10(Efit), overlays );
182  }
183 
184  }
185  }
186  }
187  STATUS("Number of events input " << setw(8) << right << job.GetBinContent(1) << endl);
188  out.Write();
189  out.Close();
190 }
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
debug
int debug
debug level
Definition: JSirene.cc:69
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition: PhysicsEvent.hh:434
JGEOMETRY3D::JVector3D::getY
double getY() const
Get y position.
Definition: JVector3D.hh:104
JGEOMETRY3D::JVector3D::getX
double getX() const
Get x position.
Definition: JVector3D.hh:94
JGEOMETRY3D::JVersor3D::getDZ
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:117
JLANG::JException
General exception.
Definition: JException.hh:24
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1714
JSUPPORT::JMultipleFileScanner::getCounter
counter_type getCounter() const
Get counter.
Definition: JMultipleFileScanner.hh:382
JTOOLS::JRange::getLowerLimit
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
JTOOLS::JRange::getUpperLimit
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
KM3NETDAQ::JDAQEventHeader::getOverlays
unsigned int getOverlays() const
Get number of overlays.
Definition: JDAQEventHeader.hh:72
KM3NETDAQ::JDAQEvent::size
unsigned int size() const
Get number of hits.
JMUONGANDALF
static const int JMUONGANDALF
Definition: definitions/reconstruction.hh:16
JMUONPREFIT
static const int JMUONPREFIT
Definition: definitions/reconstruction.hh:14
JMUONENERGY
static const int JMUONENERGY
Definition: definitions/reconstruction.hh:17
JMUONSIMPLEX
static const int JMUONSIMPLEX
Definition: definitions/reconstruction.hh:15
JMUONSTART
static const int JMUONSTART
Definition: definitions/reconstruction.hh:18
JENERGY_NUMBER_OF_HITS
static const int JENERGY_NUMBER_OF_HITS
number of hits from JEnergy.cc
Definition: fitparameters.hh:32
JGANDALF_BETA0_RAD
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v3.4.0-8-ge14cb17 https://git.km3net.de/common/km3net-dataformat.
Definition: fitparameters.hh:12
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition: JAAnetToolkit.hh:262
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
JROOT::advance
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
Definition: JROOT/JCounter.hh:35
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition: DataQueue.cc:39
std
Definition: JSTDTypes.hh:14
JACOUSTICS::JEvt
Acoustic event fit.
Definition: JAcoustics/JEvt.hh:304
JFIT::JHistory::is_event
Auxiliary class to test history.
Definition: JHistory.hh:105
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:84
JRECONSTRUCTION::JQualitySorter
General purpose sorter of fit results.
Definition: JReconstruction/JEvtToolkit.hh:288
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:45
JSUPPORT::JTriggeredFileScanner
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
Definition: JTriggeredFileScanner.hh:41
JSUPPORT::JTriggeredFileScanner::hasNext
virtual bool hasNext() override
Check availability of next element.
Definition: JTriggeredFileScanner.hh:83
JSUPPORT::JTriggeredFileScanner::next
virtual const multi_pointer_type & next() override
Get next element.
Definition: JTriggeredFileScanner.hh:103
quality_sorter
Reconstruction type dependent comparison of track quality.
Definition: tools/reconstruction.hh:85
Generated by doxygen 1.9.1