Jpp 20.0.0-27-g39925593c-D
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> quality_range;
61 JRange<double> energy_range;
62 JRange<double> coszenith_range;
63 int debug;
64
65 try {
66
67 JParser<> zap("Program to histogram fit results from reconstructed data.");
68
69 zap['f'] = make_field(inputFile);
70 zap['o'] = make_field(outputFile) = "postfit.root";
71 zap['n'] = make_field(numberOfEvents) = JLimit::max();
72 zap['N'] = make_field(numberOfPrefits) = 1;
73 zap['A'] = make_field(application) = JMUONENERGY, JMUONPREFIT, JMUONSIMPLEX, JMUONGANDALF, JMUONSTART; //makes histograms after this stage of the reconstruction
74 zap['L'] = make_field(quality_sorter) = JPARSER::initialised();
75 zap['Q'] = make_field(quality_range) = JRange<double>();
76 zap['E'] = make_field(energy_range) = JRange<double>();
77 zap['c'] = make_field(coszenith_range) = JRange<double>(-1.0,1.0);
78 zap['d'] = make_field(debug) = 2;
79
80 zap(argc, argv);
81 }
82
83 catch(const JException& error) {
84 FATAL(error);
85 }
86
87 const double rad_to_deg = 180/M_PI;
88 const double R_m = 700.0;
89 const int MAX_OVERLAYS = 300;
90 const double MAX_TRIGGERED_HITS = 2500;
91 const double E_RANGE_MIN = -1.0;
92 const double E_RANGE_MAX = 7;
93 const double Z_RANGE_MIN = coszenith_range.getLowerLimit();
94 const double Z_RANGE_MAX = coszenith_range.getUpperLimit();
95
96 TFile out(outputFile.c_str(), "recreate");
97 TH1D job("job", NULL, 100, 0.5, 100.5);
98 TH1D hz ("hz", "; cos(#theta_{zenith})" , 21 , Z_RANGE_MIN, Z_RANGE_MAX);
99 TH1D ho ("ho", "; # of overlays" , MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
100 TH2D hzo ("hzo", "; cos(#theta_{zenith}) ; # of overlays" , 21, Z_RANGE_MIN, Z_RANGE_MAX, MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
101 TH2D hxy ("hxy", "; x position ; y position" , 201, -R_m, R_m, 201, -R_m, R_m);
102 TH1D hn ("hn", "; # of triggered hits" , 600, -0.5, MAX_TRIGGERED_HITS-0.5);
103 TH1D hN ("hN", "; JEnergy # of hits" , 600, -0.5, MAX_TRIGGERED_HITS-0.5);
104 TH1D hq ("hq", "; quality parameter" , 200, -200.0, 800);
105 TH1D hb0 ("hb0", "; log(beta0)" , 60, -2.0, 3.5);
106 TH1D he ("he", "; log(Ereco [GeV])" , 75, E_RANGE_MIN, E_RANGE_MAX);
107 TH2D heo ("heo", "; log(Ereco [GeV]) ; # of overlays" , 75, E_RANGE_MIN, E_RANGE_MAX, MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
108 TH2D hen ("hen", "; log(Ereco [GeV]) ; # of triggered hits" , 75, E_RANGE_MIN, E_RANGE_MAX, 600 , -0.5, MAX_TRIGGERED_HITS-0.5);
109 TH2D heN ("heN", "; log(Ereco [GeV]) ; JEnergy # of hits" , 75, E_RANGE_MIN, E_RANGE_MAX, 600 , -0.5, MAX_TRIGGERED_HITS-0.5);
110 TH2D hzq ("hzq", "; cos(#theta_{zenith}); quality" , 21, Z_RANGE_MIN, Z_RANGE_MAX, 600, -200.0, 800.0);
111 TH2D hzn ("hzn", "; cos(#theta_{zenith}); # of triggered hits", 21, Z_RANGE_MIN, Z_RANGE_MAX, 600 , -0.5, MAX_TRIGGERED_HITS-0.5);
112 TH2D hzN ("hzN", "; cos(#theta_{zenith}); JEnergy # of hits" , 21, Z_RANGE_MIN, Z_RANGE_MAX, 600, -0.5, MAX_TRIGGERED_HITS-0.5);
113 TH2D hze ("hze", "; cos(#theta_{zenith}); log(Ereco [GeV])" , 21, Z_RANGE_MIN, Z_RANGE_MAX, 75, E_RANGE_MIN, E_RANGE_MAX);
114 TH2D hzb0("hzb0", "; cos(#theta_{zenith}); log(beta0)" , 21, Z_RANGE_MIN, Z_RANGE_MAX, 60, E_RANGE_MIN, 3.5);
115
116
117 while (inputFile.hasNext()) { // Loop over each reconstructed event
118
119 STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); // if debug > set value, print out this line
120
121 multi_pointer_type ps = inputFile.next();
122
123 JDAQEvent* tev = ps;
124 JEvt* evt = ps;
125
126 job.Fill(1.0);
127
128 if (!evt->empty()) {
129
130 JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_application(application));
131
132 if (evt->begin() == __end) {
133 continue;
134 }
135
136 if (numberOfPrefits > 0) {
137
138 JEvt::iterator __q = __end;
139
140 advance(__end = evt->begin(), min(numberOfPrefits, (size_t) distance(evt->begin(), __q)));
141
142 partial_sort(evt->begin(), __end, __q, quality_sorter);
143
144 } else {
145
146 sort(evt->begin(), __end, quality_sorter);
147 }
148
149 if (numberOfPrefits > 0) {
150 advance(__end = evt->begin(), min(numberOfPrefits, evt->size()));
151 }
152
153 for (JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) {
154
155 JTrack3D track = getTrack(*fit);
156
157 if (coszenith_range(track.getDZ()) &&
158 quality_range (fit->getQ()) &&
159 energy_range (fit->getE())) {
160
161 const double overlays = tev->getOverlays();
162 const double Efit = fit->getE();
163
164 ho .Fill(overlays);
165 hn .Fill((double) tev->size<JDAQTriggeredHit>());
166 hz .Fill(track.getDZ());
167 hzn.Fill(track.getDZ(), (double) tev->size<JDAQTriggeredHit>());
168 hzo.Fill(track.getDZ(), overlays) ;
169 hxy.Fill(track.getX(), track.getY());
170 hq .Fill(fit->getQ());
171 hzq.Fill(track.getDZ(), fit->getQ() );
172
173 if (fit->hasW(JENERGY_NUMBER_OF_HITS)) {
174 hN .Fill( fit->getW(JENERGY_NUMBER_OF_HITS) );
175 heN .Fill(log10(Efit), fit->getW(JENERGY_NUMBER_OF_HITS) );
176 hzN .Fill(track.getDZ(), fit->getW(JENERGY_NUMBER_OF_HITS) );
177 }
178
179 if (fit->hasW(JGANDALF_BETA0_RAD)) {
180 hb0 .Fill(log10( rad_to_deg* fit->getW(JGANDALF_BETA0_RAD) ));
181 hzb0 .Fill(track.getDZ(), log10( rad_to_deg* fit->getW(JGANDALF_BETA0_RAD)) );
182 }
183 he .Fill(log10(Efit));
184 hen .Fill(log10(Efit), (double) tev->size<JDAQTriggeredHit>());
185 hze .Fill(track.getDZ(), log10(Efit) );
186 heo .Fill(log10(Efit), overlays );
187 }
188
189 }
190 }
191 }
192 STATUS("Number of events input " << setw(8) << right << job.GetBinContent(1) << endl);
193 out.Write();
194 out.Close();
195}
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 build/antares-dataformat/PhysicsEvent.hh:434
JLANG::JException
General exception.
Definition JException.hh:24
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JENERGY_NUMBER_OF_HITS
static const int JENERGY_NUMBER_OF_HITS
number of hits see JRECONSTRUCTION::JMuonEnergy
Definition fitparameters.hh:33
JGANDALF_BETA0_RAD
static const int JGANDALF_BETA0_RAD
ile KM3NeT Data Definitions v3.6.2 https://git.km3net.de/common/km3net-dataformat
Definition fitparameters.hh:13
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition JAAnetToolkit.hh:278
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
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:188
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:426
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
quality_sorter
Reconstruction type dependent comparison of track quality.
Definition tools/reconstruction.hh:85
Generated by doxygen 1.10.0