Jpp  18.2.0-rc.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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

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 }
JMUONSTART
static const int JMUONSTART
Definition: definitions/reconstruction.hh:18
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1514
JLANG::JException
General exception.
Definition: JException.hh:24
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition: JAAnetToolkit.hh:259
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
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JRECONSTRUCTION::JQualitySorter
General purpose sorter of fit results.
Definition: JReconstruction/JEvtToolkit.hh:260
JSUPPORT::JTriggeredFileScanner
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
Definition: JTriggeredFileScanner.hh:39
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
KM3NETDAQ::JDAQEvent::size
unsigned int size() const
Get number of hits.
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
JMUONPREFIT
static const int JMUONPREFIT
Definition: definitions/reconstruction.hh:14
JACOUSTICS::JEvt
Acoustic event fit.
Definition: JAcoustics/JEvt.hh:299
JENERGY_NUMBER_OF_HITS
static const int JENERGY_NUMBER_OF_HITS
number of hits from JEnergy.cc
Definition: fitparameters.hh:29
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
log10
set_variable E_E log10(E_{fit}/E_{#mu})"
JFIT::JHistory::is_event
Auxiliary class to test history.
Definition: JHistory.hh:105
JMUONGANDALF
static const int JMUONGANDALF
Definition: definitions/reconstruction.hh:16
JGANDALF_BETA0_RAD
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v3.1.0-13-g917945e https://git.km3net.de/common/km3net-dataformat.
Definition: fitparameters.hh:12
KM3NETDAQ::JDAQEventHeader::getOverlays
unsigned int getOverlays() const
Get number of overlays.
Definition: JDAQEventHeader.hh:87
JGEOMETRY3D::JVector3D::getY
double getY() const
Get y position.
Definition: JVector3D.hh:104
quality_sorter
Reconstruction type dependent comparison of track quality.
Definition: tools/reconstruction.hh:85
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
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
JMUONSIMPLEX
static const int JMUONSIMPLEX
Definition: definitions/reconstruction.hh:15
JGEOMETRY3D::JVector3D::getX
double getX() const
Get x position.
Definition: JVector3D.hh:94
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:84
JGEOMETRY3D::JVersor3D::getDZ
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:117
JMUONENERGY
static const int JMUONENERGY
Definition: definitions/reconstruction.hh:17
debug
int debug
debug level
Definition: archive-put-wiki-detectors.sh:92
DEBUG
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62
Generated by   doxygen 1.8.5