Jpp  15.0.4
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JAAMuonPostfit.cc File Reference

Example program to analyse track fit results from Evt formatted data. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"
#include "km3net-dataformat/tools/reconstruction.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JPMTRouter.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JSupport/JLimit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JEvtWeightFileScannerSet.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JTools/JRange.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Example program to analyse track fit results from Evt formatted data.

Author
bofearraigh, bjung

Definition in file JAAMuonPostfit.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 40 of file JAAMuonPostfit.cc.

41 {
42  using namespace std;
43  using namespace JPP;
44  using namespace KM3NETDAQ;
45 
46  JMultipleFileScanner_t inputFiles;
47 
48  JLimit_t numberOfEvents;
49 
50  string detectorFile;
51  string outputFile;
52 
53  JRange<double> energyRange;
54  JRange<double> coszenithRange;
55  JRange<unsigned int> triggeredHitsRange;
56 
57  int weightIndex;
58  int debug;
59 
60  try {
61 
62  JProperties cuts;
63 
64  cuts.insert(gmake_property(energyRange));
65  cuts.insert(gmake_property(coszenithRange));
66  cuts.insert(gmake_property(triggeredHitsRange));
67 
68  JParser<> zap("Example program to analyse track fit results from Evt formatted data.");
69 
70  zap['f'] = make_field(inputFiles);
71  zap['n'] = make_field(numberOfEvents) = JLimit::max();
72  zap['a'] = make_field(detectorFile) = "";
73  zap['o'] = make_field(outputFile) = "histogram.root";
74  zap['@'] = make_field(cuts) = JPARSER::initialised();
75  zap['W'] = make_field(weightIndex) = -1, 2, 3, 4;
76  zap['d'] = make_field(debug) = 2;
77 
78  zap(argc, argv);
79  }
80  catch(const exception& error) {
81  FATAL(error.what() << endl);
82  }
83 
84  if (weightIndex > 0 && numberOfEvents != JLimit::max()) {
85  FATAL("Cannot apply weighting to limited number of events.");
86  }
87 
88 
89  JDetector detector;
90 
91  if (detectorFile != "") {
92  try {
93  load(detectorFile, detector);
94  }
95  catch(const JException& error) {
96  FATAL(error);
97  }
98  }
99 
100 
101  const JPMTRouter router(detector);
102 
103  // Histogram bins
104 
105  const Int_t N_Nhits = 400;
106  const Double_t Nhits_min = 0.0 - 0.5;
107  const Double_t Nhits_max = 1000.0 - 0.5;
108 
109  const Int_t N_Noverlays = 40;
110  const Double_t Noverlays_min = 0.0 - 0.5;
111  const Double_t Noverlays_max = 1000.0 - 0.5;
112 
113  const Int_t N_radius = 201;
114  const Double_t radius_min = -700.0 - 0.5;
115  const Double_t radius_max = +700.0 + 0.5;
116 
117  const Int_t N_height = 18;
118  const Double_t height_min = 20.0;
119  const Double_t height_max = 200.0;
120 
121  const Int_t N_ToTfrac = 20;
122  const Double_t ToTfrac_min = 0.0;
123  const Double_t ToTfrac_max = 1.0;
124 
125  const Int_t N_dt = 1500;
126  const Double_t dt_min = -500.0 - 0.5;
127  const Double_t dt_max = 1000.0 - 0.5;
128 
129  const Int_t N_zenith = 21;
130  const Double_t zenith_min = -1.0;
131  const Double_t zenith_max = +1.0;
132 
133  const Int_t N_energy = 60;
134  const Double_t energy_min = -2.0;
135  const Double_t energy_max = 4.0;
136 
137  const Int_t N_quality = 100;
138  const Double_t quality_min = -200.0;
139  const Double_t quality_max = +800.0;
140 
141  const Int_t N_beta0 = 60;
142  const Double_t beta0_min = -2.0;
143  const Double_t beta0_max = +3.5;
144 
145 
146  // Set of histograms for data/MC-comparisons
147 
148  TFile out(outputFile.c_str(), "recreate");
149 
150  TH1D hs ("hs", "; #snapshot hits",
151  N_Nhits, Nhits_min, Nhits_max);
152  TH1D hn ("hn", "; #triggered hits",
153  N_Nhits, Nhits_min, Nhits_max);
154  TH1D ho ("ho", "; #overlays",
155  N_Noverlays, Noverlays_min, Noverlays_max);
156 
157  TH1D hz ("hz", "; cos(#theta)",
158  N_zenith, zenith_min, zenith_max);
159  TH1D he ("he", "; log_{10}(E_{reco} / 1 GeV)",
160  N_energy, energy_min, energy_max);
161  TH1D ht ("ht", "; #Delta t [ns]",
162  N_dt, dt_min, dt_max);
163 
164  TH1D hZ ("hZ", "; longitudinal ToT-barycenter [m]",
165  N_height, height_min, height_max);
166  TH1D hT0 ("hT0", "; ToT-fraction below earliest triggered hit [ns]",
167  N_ToTfrac, ToTfrac_min, ToTfrac_max);
168  TH1D hT1 ("hT1", "; ToT-fraction above earliest triggered hit [ns]",
169  N_ToTfrac, ToTfrac_min, ToTfrac_max);
170 
171  TH1D hq ("hq", "; quality",
172  N_quality, quality_min, quality_max);
173  TH1D hb0 ("hb0", "; #beta_{0}",
174  N_beta0, beta0_min, beta0_max);
175 
176  TH2D hzo ("hzo", "; cos(#theta); #overlays",
177  N_zenith, zenith_min, zenith_max,
178  N_Noverlays, Noverlays_min, Noverlays_max);
179  TH2D hze ("hze", "; cos(#theta); log_{10}(E_{reco} / 1 GeV)",
180  N_zenith, zenith_min, zenith_max,
181  N_energy, energy_min, energy_max);
182  TH2D heo ("heo", "; log_{10}(E_{reco} / 1 GeV); #overlays",
183  N_energy, energy_min, energy_max,
184  N_Noverlays, Noverlays_min, Noverlays_max);
185  TH2D hzt ("hzt", "; cos(#theta); #Delta t [ns]",
186  N_zenith, zenith_min, zenith_max,
187  N_dt, dt_min, dt_max);
188  TH2D het ("het", "; log_{10}(E_{reco} / 1 GeV); #Delta t [ns]",
189  N_energy, energy_min, energy_max,
190  N_dt, dt_min, dt_max);
191  TH2D hot ("hot", "; #overlays; #Delta t [ns]",
192  N_Noverlays, Noverlays_min, Noverlays_max,
193  N_dt, dt_min, dt_max);
194 
195  TH2D hzq ("hzq", "; cos(#theta); quality",
196  N_zenith, zenith_min, zenith_max,
197  N_quality, quality_min, quality_max);
198  TH2D heq ("heq", "; log_{10}(E_{reco} / 1 GeV); quality",
199  N_energy, energy_min, energy_max,
200  N_quality, quality_min, quality_max);
201  TH2D hoq ("hoq", "; #overlays; quality",
202  N_Noverlays, Noverlays_min, Noverlays_max,
203  N_quality, quality_min, quality_max);
204 
205  TH2D hob0("hob0", "; #overlays; log_{10}(#beta_{0})",
206  N_Noverlays, Noverlays_min, Noverlays_max,
207  N_beta0, beta0_min, beta0_max);
208  TH2D heb0("heb0", "; log_{10}(E_{reco} / 1 GeV); log_{10}(#beta_{0})",
209  N_energy, energy_min, energy_max,
210  N_beta0, beta0_min, beta0_max);
211  TH2D hzb0("hzb0", "; cos(#theta); log_{10}(#beta_{0})",
212  N_zenith, zenith_min, zenith_max,
213  N_beta0, beta0_min, beta0_max);
214 
215  TH2D hxy ("hxy", "; x [m]; y [m]",
216  N_radius, radius_min, radius_max,
217  N_radius, radius_min, radius_max);
218 
219 
220  JEvtWeightFileScannerSet<> scanners(inputFiles);
221 
222  STATUS(RIGHT(15) << "header ID" << RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
223 
224  for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
225 
226  scanner->setLimit(numberOfEvents);
227 
228  while (scanner->hasNext()) {
229 
230  const Evt* evt = scanner->next();
231 
232  if (has_reconstructed_jppmuon(*evt)) {
233 
234  DEBUG(get_best_reconstructed_jppmuon(*evt).comment << endl);
235 
236  const Trk& best_trk = get_best_reconstructed_jppmuon(*evt);
237  const JTrack3E& tb = getTrack(best_trk);
238 
239  const double logE = log10(tb.getE());
240  const size_t NsnapshotHits = evt->hits.size();
241  const size_t NtriggeredHits = count_if(evt->hits.begin(), evt->hits.end(),
242  make_predicate(&Hit::trig, (long long unsigned int) 0, JComparison::gt()));
243 
244  if (triggeredHitsRange(NtriggeredHits) &&
245  coszenithRange (tb.getDZ()) &&
246  energyRange (tb.getE())) { // Apply cuts
247 
248  // Extract event-weight
249 
250  double weight;
251 
252  if (weightIndex < 2) {
253 
254  weight = 1.0;
255 
256  } else if (weightIndex < 3) {
257 
258  weight = scanner->getWeight(*evt);
259 
260  } else if (weightIndex <= (int) evt->w.size()) {
261 
262  weight = evt->w.at(weightIndex - 1);
263 
264  } else {
265 
266  WARNING('w' << weightIndex << "-weight does not exist for event #" << scanner->getCounter() << "; skip." << endl);
267  continue;
268  }
269 
270  STATUS(RIGHT (15) << distance(scanners.begin(), scanner) <<
271  RIGHT (15) << scanner->getCounter() <<
272  SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);
273 
274  // Fill histograms
275 
276  hs .Fill(NsnapshotHits, weight);
277  hn .Fill(NtriggeredHits, weight);
278  ho .Fill(evt->overlays, weight);
279 
280  hz .Fill(tb.getDZ(), weight);
281  he .Fill(logE, weight);
282  hq .Fill(best_trk.lik, weight);
283 
284  hzo.Fill(tb.getDZ(), evt->overlays, weight);
285  hze.Fill(tb.getDZ(), logE, weight);
286  heo.Fill(logE, evt->overlays, weight);
287 
288  hoq.Fill(evt->overlays, best_trk.lik, weight);
289  hzq.Fill(tb.getDZ(), best_trk.lik, weight);
290  heq.Fill(logE, best_trk.lik, weight);
291 
292  hxy.Fill(tb.getX(), tb.getY(), weight);
293 
294  if (best_trk.fitinf[JGANDALF_BETA0_RAD]) {
295 
296  const double logb0 = log10(180.0 / M_PI * best_trk.fitinf[JGANDALF_BETA0_RAD]);
297 
298  hb0 .Fill(logb0, weight);
299  hob0.Fill(evt->overlays, logb0, weight);
300  heb0.Fill(logE, logb0, weight);
301  hzb0.Fill(tb.getDZ(), logb0, weight);
302  }
303 
304  double ToTbelow = 0.0;
305  double ToTabove = 0.0;
306  double ToTtotal = 0.0;
307  double ToTweightedZ = 0.0;
308 
309  vector<Hit>::const_iterator firstHit = min_element(evt->hits.cbegin(), evt->hits.cend(),
310  make_comparator(&Hit::t));
311 
312  for (vector<Hit>::const_iterator hit = evt->hits.cbegin(); hit != evt->hits.cend(); ++hit) {
313 
314  if (router.hasPMT(hit->pmt_id)) {
315 
316  const JPMT& pmt = router.getPMT(hit->pmt_id);
317 
318  const double t0 = tb.getT(pmt);
319  const double t1 = getTime(*hit);
320 
321  ht .Fill(t1 - t0, weight);
322  hot.Fill(evt->overlays, t1 - t0, weight);
323  hzt.Fill(tb.getDZ(), t1 - t0, weight);
324  het.Fill(log10(best_trk.E), t1 - t0, weight);
325 
326  if (hit->trig > 0) {
327 
328  ToTtotal += hit->tot;
329  ToTweightedZ += hit->tot * hit->pos.z;
330 
331  if (hit->pos.z < firstHit->pos.z) {
332  ToTbelow += hit->tot;
333  } else if (hit->pos.z >= firstHit->pos.z) {
334  ToTabove += hit->tot;
335  }
336  }
337  }
338  }
339 
340  hT0.Fill(ToTbelow / ToTtotal, weight);
341  hT1.Fill(ToTabove / ToTtotal, weight);
342  hZ .Fill(ToTweightedZ / ToTtotal, weight);
343  }
344  }
345  }
346  }
347 
348  STATUS(endl);
349 
350  out.Write();
351  out.Close();
352 }
JDETECTOR::JPMTRouter
Router for direct addressing of PMT data in detector data structure.
Definition: JPMTRouter.hh:33
JPARSER::JParser
Utility class to parse command line options.
Definition: JParser.hh:1500
JLANG::JException
General exception.
Definition: JException.hh:23
WARNING
#define WARNING(A)
Definition: JMessage.hh:65
JLANG::make_predicate
JPredicate< JResult_t T::*, JComparison::eq > make_predicate(JResult_t T::*member, const JResult_t value)
Helper method to create predicate for data member.
Definition: JPredicate.hh:128
JLANG::JComparison::gt
Greater than.
Definition: JComparison.hh:73
JLANG::make_comparator
JComparator< JResult_t T::*, JComparison::lt > make_comparator(JResult_t T::*member)
Helper method to create comparator between values of data member.
Definition: JComparator.hh:185
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition: JAAnetToolkit.hh:246
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
gmake_property
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
Definition: JProperties.hh:1291
STATUS
#define STATUS(A)
Definition: JMessage.hh:63
JDETECTOR::JDetector
Detector data structure.
Definition: JDetector.hh:89
Trk::comment
std::string comment
use as you like
Definition: Trk.hh:34
JEEP::JProperties
Utility class to parse parameter values.
Definition: JProperties.hh:496
Evt::w
std::vector< double > w
MC: Weights w[0]=w1, w[1]=w2, w[2]]=w3 (see e.g. Tag list)
Definition: Evt.hh:39
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
JAANET::getTime
double getTime(const Hit &hit)
Get true time of hit.
Definition: JAAnetToolkit.hh:84
outputFile
string outputFile
Definition: JDAQTimesliceSelector.cc:37
has_reconstructed_jppmuon
bool has_reconstructed_jppmuon(const Evt &evt)
Test whether given event has a track with muon reconstruction.
Definition: tools/reconstruction.hh:305
Trk::E
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
JGEOMETRY3D::JTrack3E
3D track with energy.
Definition: JTrack3E.hh:30
JSUPPORT::JEvtWeightFileScannerSet::iterator
std::vector< value_type >::iterator iterator
Definition: JEvtWeightFileScannerSet.hh:57
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
Evt::overlays
unsigned int overlays
number of overlaying triggered events
Definition: Evt.hh:31
JAANET::detector
Detector file.
Definition: JHead.hh:196
JGEOMETRY3D::JTrack3E::getE
double getE() const
Get energy.
Definition: JTrack3E.hh:93
get_best_reconstructed_jppmuon
const Trk & get_best_reconstructed_jppmuon(const Evt &evt)
Get best reconstructed muon.
Definition: tools/reconstruction.hh:387
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
JDETECTOR::JPMT
Data structure for PMT geometry, calibration and status.
Definition: JPMT.hh:43
Trk::fitinf
std::vector< double > fitinf
place to store additional fit info, for jgandalf, see JFitParameters.hh
Definition: Trk.hh:31
JGANDALF_BETA0_RAD
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v2.1.0-7-g97845ea https://git.km3net.de/common/km3net-dataformat.
Definition: fitparameters.hh:12
JGEOMETRY3D::JVector3D::getY
double getY() const
Get y position.
Definition: JVector3D.hh:104
debug
int debug
debug level
Definition: JSirene.cc:63
JGEOMETRY3D::JTrack3D::getT
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JTrack3D.hh:126
FATAL
#define FATAL(A)
Definition: JMessage.hh:67
JSUPPORT::JMultipleFileScanner_t
Auxiliary base class for list of file names.
Definition: JMultipleFileScanner.hh:45
JDETECTOR::load
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Definition: JDetectorToolkit.hh:475
Hit::trig
ULong64_t trig
non-zero if the hit is a trigger hit.
Definition: Hit.hh:18
Trk::lik
double lik
likelihood or lambda value (for aafit, lambda)
Definition: Trk.hh:23
Hit::t
double t
hit time (from tdc+calibration or MC truth)
Definition: Hit.hh:23
DEBUG
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
Definition: JCalibratePMT.sh:21
JGEOMETRY3D::JVector3D::getX
double getX() const
Get x position.
Definition: JVector3D.hh:94
JSUPPORT::JEvtWeightFileScannerSet
Auxiliary class for organising Monte Carlo file scanners associated with event weighters.
Definition: JEvtWeightFileScannerSet.hh:51
Evt::hits
std::vector< Hit > hits
list of hits
Definition: Evt.hh:35
detector
do set_variable DETECTOR_TXT $WORKDIR detector
Definition: detector-XY:fit1d.sh:38
Trk
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:484
JGEOMETRY3D::JVersor3D::getDZ
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:117
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19
JTRIGGER::weight
std::vector< double > weight
Definition: JAlgorithm.hh:417
Generated by   doxygen 1.8.5