Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
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_application(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:72
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
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:1698
JSUPPORT::JMultipleFileScanner::getCounter
counter_type getCounter() const
Get counter.
Definition JMultipleFileScanner.hh:382
JTOOLS::JRange
Range of values.
Definition JRange.hh:42
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.5.1-1-gd514d72 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
JACOUSTICS::JEvt
Acoustic event fit.
Definition JAcoustics/JEvt.hh:310
JFIT::JHistory::is_application
Auxiliary class to test history.
Definition JHistory.hh:167
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
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::JLimit::max
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
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.12.0