Jpp  15.0.4
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JAAMuonPostfit.cc
Go to the documentation of this file.
1 
2 #include <string>
3 #include <iostream>
4 #include <iomanip>
5 #include <cmath>
6 
7 #include "TROOT.h"
8 #include "TFile.h"
9 #include "TH1D.h"
10 
15 
16 #include "Jeep/JParser.hh"
17 #include "Jeep/JMessage.hh"
18 
19 #include "JDetector/JDetector.hh"
20 #include "JDetector/JPMTRouter.hh"
22 
23 #include "JAAnet/JHeadToolkit.hh"
24 #include "JAAnet/JAAnetToolkit.hh"
25 
26 #include "JSupport/JLimit.hh"
27 #include "JSupport/JSupport.hh"
31 
32 #include "JTools/JRange.hh"
33 
34 /**
35  * \file
36  *
37  * Example program to analyse track fit results from Evt formatted data.
38  * \author bofearraigh, bjung
39  */
40 int main(int argc, char **argv)
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 
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 
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(),
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 }
Router for direct addressing of PMT data in detector data structure.
Definition: JPMTRouter.hh:33
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
#define WARNING(A)
Definition: JMessage.hh:65
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
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Greater than.
Definition: JComparison.hh:73
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
JTrack3E getTrack(const Trk &track)
Get track.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
#define gmake_property(A)
macro to convert (template) parameter to JPropertiesElement object
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:89
Auxiliary methods for selection of reconstructed tracks.
std::string comment
use as you like
Definition: Trk.hh:34
Utility class to parse parameter values.
Definition: JProperties.hh:496
bool hasPMT(const JObjectID &id) const
Has PMT.
Definition: JPMTRouter.hh:114
const JPMT & getPMT(const JPMTAddress &address) const
Get PMT.
Definition: JPMTRouter.hh:90
std::vector< double > w
MC: Weights w[0]=w1, w[1]=w2, w[2]]=w3 (see e.g. Tag list)
Definition: Evt.hh:39
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:66
double getTime(const Hit &hit)
Get true time of hit.
string outputFile
bool has_reconstructed_jppmuon(const Evt &evt)
Test whether given event has a track with muon reconstruction.
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
3D track with energy.
Definition: JTrack3E.hh:30
Data structure for detector geometry and calibration.
std::vector< value_type >::iterator iterator
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
unsigned int overlays
number of overlaying triggered events
Definition: Evt.hh:31
Detector file.
Definition: JHead.hh:196
double getE() const
Get energy.
Definition: JTrack3E.hh:93
const Trk & get_best_reconstructed_jppmuon(const Evt &evt)
Get best reconstructed muon.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1961
Data structure for PMT geometry, calibration and status.
Definition: JPMT.hh:43
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
std::vector< double > fitinf
place to store additional fit info, for jgandalf, see JFitParameters.hh
Definition: Trk.hh:31
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v2.1.0-7-g97845ea https://git.km3net.de/common/km3net-dataformat.
double getY() const
Get y position.
Definition: JVector3D.hh:104
Direct access to PMT in detector data structure.
int debug
debug level
Definition: JSirene.cc:63
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JTrack3D.hh:126
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Auxiliary base class for list of file names.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
ULong64_t trig
non-zero if the hit is a trigger hit.
Definition: Hit.hh:18
double lik
likelihood or lambda value (for aafit, lambda)
Definition: Trk.hh:23
Utility class to parse command line options.
double t
hit time (from tdc+calibration or MC truth)
Definition: Hit.hh:23
then usage $script< input_file >< detector_file > fi set_variable OUTPUT_DIR set_variable SELECTOR JDAQTimesliceL1 set_variable DEBUG case set_variable DEBUG
double getX() const
Get x position.
Definition: JVector3D.hh:94
Auxiliary class for organising Monte Carlo file scanners associated with event weighters.
std::vector< Hit > hits
list of hits
Definition: Evt.hh:35
do set_variable DETECTOR_TXT $WORKDIR detector
Auxiliaries for defining the range of iterations of objects.
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:484
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:117
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:19
std::vector< double > weight
Definition: JAlgorithm.hh:417