Jpp test-rotations-new
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 JMUONGANDALF
static const int JMUONPREFIT
static const int JMUONENERGY
static const int JMUONSIMPLEX
static const int JMUONSTART
static const int JENERGY_NUMBER_OF_HITS
number of hits from JEnergy.cc
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v3.6.0 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:167
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.