Jpp 20.0.0-rc.2
the software that should make you happy
Loading...
Searching...
No Matches
JDataPostfit.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4
5#include "TROOT.h"
6#include "TFile.h"
7#include "TH1D.h"
8#include "TH2D.h"
9#include "TMath.h"
10
16
19
20#include "JAAnet/JHead.hh"
22#include "JAAnet/JVolume.hh"
23
27
30#include "JSupport/JSupport.hh"
31
32#include "JDAQ/JDAQEventIO.hh"
33
34#include "Jeep/JParser.hh"
35#include "Jeep/JMessage.hh"
36
37/**
38 * \file
39 *
40 * Program to histogram fit results from reconstructed data
41 * \author bofearraigh, rgracia
42 */
43int main(int argc, char **argv)
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
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
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}
string outputFile
int main(int argc, char **argv)
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
ROOT TTree parameter settings of various packages.
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
double getY() const
Get y position.
Definition JVector3D.hh:104
double getX() const
Get x position.
Definition JVector3D.hh:94
double getDZ() const
Get z direction.
Definition JVersor3D.hh:117
General exception.
Definition JException.hh:24
Utility class to parse command line options.
Definition JParser.hh:1698
counter_type getCounter() const
Get counter.
Range of values.
Definition JRange.hh:42
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
unsigned int getOverlays() const
Get number of overlays.
unsigned int size() const
Get number of hits.
static const int JENERGY_NUMBER_OF_HITS
number of hits from JMuonEnergy
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v3.6.1-2-g905a24d https://git.km3net.de/common/km3net-dataformat.
JTrack3E getTrack(const Trk &track)
Get track.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
counter_type advance(counter_type &counter, const counter_type value, const counter_type limit=std::numeric_limits< counter_type >::max())
Advance counter.
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
Acoustic event fit.
Auxiliary class to test history.
Definition JHistory.hh:188
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
General purpose sorter of fit results.
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
virtual bool hasNext() override
Check availability of next element.
virtual const multi_pointer_type & next() override
Get next element.
Reconstruction type dependent comparison of track quality.