Jpp master_rocky-44-g75b7c4f75
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JShowerPostfit.cc File Reference

Example program to histogram shower fit results: it handles both neutrino and muon productions. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile.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 "km3net-dataformat/tools/time_converter.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JVolume.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JPhysics/JConstants.hh"
#include "JTools/JQuantile.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JLang/JPredicate.hh"
#include "JPhysics/JGeanz.hh"
#include "JGeometry3D/JShower3E.hh"
#include "JGeometry3D/JCylinder3D.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

Example program to histogram shower fit results: it handles both neutrino and muon productions.

Author
mdejong, adomi

Definition in file JShowerPostfit.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 79 of file JShowerPostfit.cc.

80{
81 using namespace std;
82 using namespace JPP;
83 using namespace KM3NETDAQ;
84
85 typedef JTriggeredFileScanner<JEvt> JTriggeredFileScanner_t;
86 typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
87
88 JTriggeredFileScanner_t inputFile;
89 JLimit_t& numberOfEvents = inputFile.getLimit();
90 string outputFile;
91 size_t numberOfPrefits;
92 JQualitySorter quality_sorter;
93 JAtmosphericMuon atmosphere;
94 double Emin_GeV;
95 double NPE;
96 int application;
97 string option;
98 bool isMuon;
99 bool wrtNeutrino;
100 double radius;
101 JRange<double> height;
102 int debug;
103
104 try {
105
106 JParser<> zap("Example program to histogram fit results.");
107
108 zap['f'] = make_field(inputFile);
109 zap['o'] = make_field(outputFile) = "postfit.root";
110 zap['n'] = make_field(numberOfEvents) = JLimit::max();
111 zap['N'] = make_field(numberOfPrefits) = 1;
112 zap['L'] = make_field(quality_sorter) = JPARSER::initialised();
113 zap['E'] = make_field(Emin_GeV) = 0.0;
114 zap['M'] = make_field(NPE) = 0.0;
115 zap['A'] = make_field(application) = JSHOWERPREFIT, JSHOWERPOINTSIMPLEX, JSHOWERPOSITIONFIT, JSHOWERDIRECTIONPREFIT, JSHOWERCOMPLETEFIT, JSHOWER_BJORKEN_Y;
116 zap['a'] = make_field(atmosphere) = JAtmosphericMuon(90.0, 90.0);
117 zap['O'] = make_field(option) = "E", "N", "LINE", "LOGE";
118 zap['I'] = make_field(isMuon);
119 zap['w'] = make_field(wrtNeutrino);
120 zap['r'] = make_field(radius) = numeric_limits<double>::max(); // no containment
121 zap['z'] = make_field(height) = JRange<double>(); // no containment
122 zap['d'] = make_field(debug) = 2;
123
124 zap(argc, argv);
125 }
126 catch(const exception& error) {
127 FATAL(error.what() << endl);
128 }
129
130 cout << "APPLICATION " << application << endl;
131
132 JHead head;
133
134 try {
135 head = getCommonHeader(inputFile);
136 }
137 catch(const exception& error) {
138 FATAL(error.what() << endl);
139 }
140
141 JVolume volume(head, option != "LINE");
142 JPosition3D offset = getPosition(getOffset(head));
143 JPosition3D origin = getPosition(getOrigin(head));
144 JCylinder3D cylinder = getCylinder(head);
145
146 cylinder.add(offset);
147
148 JCylinder3D containment(JCircle2D(JPosition2D(), radius), height.getLowerLimit(), height.getUpperLimit());
149
150 containment.add(origin);
151
152 NOTICE("Offset: " << offset << endl);
153 NOTICE("Origin: " << origin << endl);
154 NOTICE("Cylinder: " << cylinder << endl);
155 NOTICE("Containment: " << containment << endl);
156
157 const double EMIN_GEV = 10e3; // MUPAGE
158
159 TFile out(outputFile.c_str(), "recreate");
160
161 TH1D job("job", NULL, 100, 0.5, 100.5);
162
163 TH1D hn("hn", "NDF (Number of Degrees of Freedom)", 250, -0.5, 249.5);
164 TH1D hq("hq", "Fit Quality", 300, 0.0, 600.0);
165 TH1D hi("hi", "Fit Index vs Best Resolution", 101, -0.5, 100.5); // to see if the fit with the best quality has also the best resolution
166 TH1D hd("hd", "Distance between Reco and True position", 2000, 0.0, 100.0); // [m]
167 TH1D hz("hz", "Longitudinal Distance between Reco and MC lepton position", 100, -200.0, 200.0); // [m]
168 TH1D ht("ht", "Time difference between Reco and MC lepton", 100, 0, 100.0); // [ns]
169 TH1D h4D("h4D", "4D-distance between reco and true vertex", 2000, 0.0, 100.0); //[m]
170
171 TH1D ha("ha", "Angle between reco direction and true direction", 1000, 0.0, 180.0); // [ns]
172
173 TH1D e0("e0", "True lepton energy", 100, volume.getXmin(), volume.getXmax()); // [log(E)]
174 TH1D e2("e2", "Reconstructed energy", 100, volume.getXmin(), volume.getXmax()); // [log(E)]
175 TH1D er("er", "Ratio of reconstructed energy and true energy", 100, -5.0, +5.0); // [log(E/E)]
176 TH1D ea("ea", "Distribution of Energy error for all the events; E_{reco}-E_{true}", 100, -30, +30);
177 TH1D ed3_5GeV("ed3_5GeV", "Distribution of Energy error for events with MC energy #in [3,5] GeV; E_{reco}-E_{true}", 100, -30, +30);
178 TH1D ed8_11GeV("ed8_11GeV", "Distribution of Energy error for events with MC energy #in [8,11] GeV; E_{reco}-E_{true}", 100, -30, +30);
179 TH1D ed15_21GeV("ed15_21GeV", "Distribution of Energy error for events with MC energy #in [15,21] GeV; E_{reco}-E_{true}", 100, -30, +30);
180 TH2D ee("ee", "; E_{true} [GeV]; E_{reco} [GeV]",
181 40, 0, 4,
182 40, 0, 4);
183 TH1D eeres("eeres","Median reco energy as a function of true energy",40,0,4);
184 TH1D eeres16("eeres16","16%Quantile reco energy as a function of true energy",40,0,4);
185 TH1D eeres84("eeres84","84%Quantile reco energy as a function of true energy",40,0,4);
186 const Int_t ny = 28800;
187 const Double_t ymin = 0.0; // [deg]
188 const Double_t ymax = 180.0; // [deg]
189
190 vector<double> X;
191
192 if (option.find('E') != string::npos) {
193
194 for (double x = volume.getXmin(); x <= volume.getXmax(); x += (volume.getXmax() - volume.getXmin()) / 20) {
195 X.push_back(x);
196 }
197
198 } else {
199
200 double x = -0.5;
201
202 for ( ; x <= 15.5; x += 1.0) { X.push_back(x); }
203 for ( ; x <= 25.5; x += 2.0) { X.push_back(x); }
204 for ( ; x <= 50.5; x += 5.0) { X.push_back(x); }
205 for ( ; x <= 100.5; x += 10.0) { X.push_back(x); }
206 for ( ; x <= 250.5; x += 20.0) { X.push_back(x); }
207 }
208
209 TH2D h2("h2", NULL, X.size() - 1, X.data(), ny, ymin, ymax);
210
211 TProfile herrorE("herrorE", "Energy Error", X.size() - 1, X.data());
212 TProfile hfracE( "hfracE", "Fractional Energy Error", X.size() - 1, X.data());
213 TProfile he("he","Reconstruction Efficiency", X.size() - 1, X.data());
214 TProfile htheta_nu_lep("htheta_nu_lep", "Angle between neutrino and primary lepton", X.size() - 1, X.data());
215 TH2D hgetln("hgetln","Kinematic angle", 40,0,4,1000,0,180);
216 TH1D hln1d("hln1d","Angle between neutrino and lepton",40,0,4);
217 TH2D hb("hb", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
218 TH2S hmca("hmca", NULL, X.size() - 1, X.data(), 1000, 0, 180);
219 TH2D hae("hae",NULL,40,0,4,1000,0,180);
220 TH1D hmae("hmae","Mean angle deviation as function of log MC energy",40,0,4);
221 TH2D hzenith("hzenith", "Reco Zenith vs MC Energy", 100, 0, 100, ny, ymin, ymax);
222 TH2D hY("hY", "Reco Bjorken Y vs MC Bjorken Y", 40, 0, 1, 40, 0, 1);
223 TH2D hby3_5GeV("hby3_5GeV", "Reco Bjorken Y vs MC Bjorken Y for events with Reco E #in [3, 5] GeV", 40, 0, 1, 40, 0, 1);
224 TH2D hby8_10GeV("hby8_10GeV", "Reco Bjorken Y vs MC Bjorken Y for events with Reco E #in [8, 10] GeV", 40, 0, 1, 40, 0, 1);
225
226 JQuantile Q("Angle", true);
227 JQuantile O("Omega", true);
228 JQuantile QE("Energy", true);
229 JQuantile Qp("Pos",true);
230 JQuantile Qt("Time",true);
231 JQuantile Q4d("4D",true);
232 while (inputFile.hasNext()) {
233
234 STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
235
236 multi_pointer_type ps = inputFile.next();
237
238 JDAQEvent* tev = ps;
239 JEvt* evt = ps; // Reco event
240 Evt* event = ps; // MC event
241
242 const time_converter converter(*event, *tev);
243
244 job.Fill(1.0);
245
246 double Enu = 0.0;
247 double Elep = 0.0;
248 double Emax = 0.0;
249
250 Trk neutrino;
251
252 if(has_neutrino(*event)){
253
254 neutrino = get_neutrino(*event);
255 Enu = neutrino.E;
256
257 }
258
259 vector<Trk>::const_iterator lepton = event->mc_trks.end(); // can be electron or muon
260
261 for (vector<Trk>::const_iterator mc_evt = event->mc_trks.begin(); mc_evt != event->mc_trks.end(); ++mc_evt) {
262
263 if (!isMuon && is_electron(*mc_evt)) {
264
265 Elep += mc_evt->E;
266
267 if (mc_evt->E > Emax) {
268 lepton = mc_evt;
269 Emax = mc_evt->E;
270 }
271
272
273 } else if(isMuon && is_muon(*mc_evt)){
274
275 Elep += mc_evt->E;
276 if (mc_evt->E > Emax) {
277 lepton = mc_evt;
278 Emax = mc_evt->E;
279 }
280 }
281 }
282
283 if (lepton == event->mc_trks.end()) {
284 continue;
285 }
286
287 job.Fill(3.0);
288
289 double true_BjY = (Enu - Elep) / Enu;
290
291 // abscissa
292
293 Double_t x = 0.0;
294
295 if (option.find('E') != string::npos){
296
297 if(wrtNeutrino == true && !isMuon) x = volume.getX(neutrino.E);
298 else if(!wrtNeutrino && !isMuon) x = volume.getX(lepton->E);
299
300 } else{
301 x = getCount(tev->begin<JDAQTriggeredHit>(), tev->end<JDAQTriggeredHit>());
302 }
303
304 // weight for efficiency determination
305 Double_t W = 0.0;
306
307 if (!evt->empty()) {
308
309 JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_event(application));
310
311 if (evt->begin() == __end) {
312 continue;
313 }
314
315 job.Fill(4.0);
316
317 //if (numberOfPrefits > 0 && numberOfPrefits < (size_t) distance(evt->begin(), evt->end())) {
318 if (numberOfPrefits > 0 ) {
319
320 JEvt::iterator __q = __end;
321
322 advance(__end = evt->begin(), min(numberOfPrefits, (size_t) distance(evt->begin(), __q)));
323
324 partial_sort(evt->begin(), __end, __q, quality_sorter);
325
326 } else {
327
328 sort(evt->begin(), __end, quality_sorter);
329 }
330
331 JEvt::iterator best = evt->begin(); // the best shower is the one with the best quality
332
333 JPosition position(getPosition(lepton->pos)); // by default the lepton pos is the interaction vertex
334 JPointing pointing(getDirection(*lepton));
335 JShowerEnergy energy(lepton->E);
336
337 JVector3D showerBrightPoint(getPosition(lepton->pos) + offset + getPosition(lepton->dir) * geanz.getMaximum(lepton->E));
338
339 if(!wrtNeutrino && !isMuon){
340 position = JPosition(showerBrightPoint); // wrt the em shower brightest point
341 } else if(wrtNeutrino == true && !isMuon){
342 position = JPosition(getPosition(neutrino));
343 pointing = JPointing(getDirection(neutrino));
344 energy = JShowerEnergy(Enu);
345 }
346
347 JEvt::iterator fit_with_smallest_error = best;
348
349 if(application == JSHOWERPREFIT || application == JSHOWERPOINTSIMPLEX || application == JSHOWERPOSITIONFIT){
350 fit_with_smallest_error = position(evt->begin(), __end);
351 } else if (application == JSHOWERENERGYPREFIT){
352 fit_with_smallest_error = energy(evt->begin(), __end);
353 } else {
354 fit_with_smallest_error = pointing(evt->begin(), __end);
355 }
356 best = fit_with_smallest_error;
357
358 const Double_t beta = pointing.getAngle(*best);
359 const double Efit = best->getE();
360
361 if (containment.is_inside(getPosition(*best))){
362
363 // selection of fit result
364 bool ok = (Efit >= Emin_GeV);
365
366 if (ok) {
367
368 W = 1.0;
369
370 job.Fill(5.0);
371
372 hn.Fill((Double_t) best->getNDF());
373 hq.Fill(best->getQ());
374 hi.Fill((Double_t) distance(evt->begin(), fit_with_smallest_error));
375 Q.put(beta);
376
377 JShower3E ta;
378
379 if(wrtNeutrino){
380 ta = getTrack(neutrino);
381 } else {
382 ta = getTrack(*lepton); // MC event
383 }
384
385 JShower3E tb = getTrack(*best); // Reco event
386
387 double zenith = tb.getDirection().getTheta() * 180 / PI;
388
389 // take into account the difefrent MC/Reco geometry and time
390 ta.add(offset);
391 tb.sub(converter.putTime());
392 double time_true = ta.getT();
393 double time_reco=tb.getT();
394 double distance_m;
395
396 if(!isMuon && !wrtNeutrino){
397 JPosition3D posDif = tb.getPosition() - JPosition3D(showerBrightPoint);
398 distance_m = posDif.getLength();
399 time_true += geanz.getMaximum(lepton->E)*getInverseSpeedOfLight();
400 hd.Fill(fabs(distance_m));
401 ht.Fill(fabs(time_reco - time_true));
402 Qp.put(fabs(distance_m));
403 Qt.put(fabs(time_reco-time_true));
404 double distance4d = sqrt(posDif.getLengthSquared() + pow(getSpeedOfLight()/getIndexOfRefraction()*(time_true-time_reco),2));
405 h4D.Fill(distance4d);
406 Q4d.put(distance4d);
407 } else {
408 JPosition3D posDif = tb.getPosition() - ta.getPosition();
409 distance_m = posDif.getLength();
410 hd.Fill(fabs(distance_m));
411 ht.Fill(fabs(time_reco - time_true));
412 Qp.put(fabs(distance_m));
413 Qt.put(fabs(time_reco-time_true));
414 double distance4d = sqrt(posDif.getLengthSquared() + pow(getSpeedOfLight()/getIndexOfRefraction()*(time_true-time_reco),2));
415 h4D.Fill(distance4d);
416 Q4d.put(distance4d);
417 }
418
419 if (has_neutrino(*event)) {
420
421 job.Fill(6.0);
422
423 JPosition3D mc_vx = getPosition(neutrino); // mc vertex
424
425 mc_vx.add(offset);
426
427 if (cylinder.is_inside(mc_vx)) {
428
429 job.Fill(7.0);
430
431 hz.Fill(tb.getIntersection(mc_vx));
432 }
433 }
434
435 if (best->getE() >= EMIN_GEV) {
436 hb.Fill(JVector2D(best->getX() - origin.getX(), best->getY() - origin.getY()).getLengthSquared(), best->getZ());
437 }
438
439 ha.Fill(getAngle(ta.getDirection(), tb.getDirection()));
440
441 htheta_nu_lep.Fill(x, getAngle(getDirection(neutrino), getDirection(*lepton)));
442 hgetln.Fill(log10(Enu), getAngle(getDirection(neutrino), getDirection(*lepton)));
443 for (int i = 1; i <= hgetln.GetNbinsX(); ++i) {
444 std::vector<double> deltaThetaValues;
445 for (int j = 1; j <= hgetln.GetNbinsY(); ++j) {
446 double binContent = hgetln.GetBinContent(i, j);
447 // Add Delta theta values based on the bin content (number of events)
448 for (int k = 0; k < binContent; ++k) {
449 double deltaTheta = hgetln.GetYaxis()->GetBinCenter(j);
450 deltaThetaValues.push_back(deltaTheta);
451 }
452 }
453
454 // Calculate the median using TMath::Median
455 double medianDeltaTheta = TMath::Median(deltaThetaValues.size(), &deltaThetaValues[0]);
456
457 // Set the median Delta theta as a function of energy in the graph
458 hln1d.SetBinContent(i, medianDeltaTheta);
459 }
460 hmca.Fill(x, getAngle(ta.getDirection(), tb.getDirection()));
461 hae.Fill(log10(Enu),getAngle(ta.getDirection(), tb.getDirection()));
462 for (int i = 1; i <= hae.GetNbinsX(); ++i) {
463 std::vector<double> deltaThetaValues;
464 for (int j = 1; j <= hae.GetNbinsY(); ++j) {
465 double binContent = hae.GetBinContent(i, j);
466 // Add Delta theta values based on the bin content (number of events)
467 for (int k = 0; k < binContent; ++k) {
468 double deltaTheta = hae.GetYaxis()->GetBinCenter(j);
469 deltaThetaValues.push_back(deltaTheta);
470 }
471 }
472
473 // Calculate the median using TMath::Median
474 double medianDeltaTheta = TMath::Median(deltaThetaValues.size(), &deltaThetaValues[0]);
475
476 // Set the median Delta theta as a function of energy in the graph
477 hmae.SetBinContent(i, medianDeltaTheta);
478 }
479 if(application == JSHOWER_BJORKEN_Y) {
480
481 hY.Fill(true_BjY, best->getW(5));
482
483 if(Efit > 2.5 && Efit < 6) hby3_5GeV.Fill(true_BjY, best->getW(5));
484 else if(Efit > 7.5 && Efit < 11) hby8_10GeV.Fill(true_BjY, best->getW(5));
485
486 }
487 if(wrtNeutrino){
488
489 e0.Fill(volume.getX(Enu, true));
490 er.Fill(volume.getX(Efit) - volume.getX(Enu));
491 ee.Fill(volume.getX(Enu), volume.getX(Efit));
492 ea.Fill(Efit - Enu);
493 hzenith.Fill(Enu, zenith);
494 QE.put(log10(Efit/Enu));
495 herrorE.Fill(x, volume.getX(Efit) - volume.getX(Enu));
496 hfracE.Fill(x, fabs(volume.getX(Efit) - volume.getX(Enu))/volume.getX(Enu));
497
498 if(Enu >= 3 && Enu <= 5) ed3_5GeV.Fill(Efit - Enu);
499 else if(Enu >= 8 && Enu <= 11) ed8_11GeV.Fill(Efit - Enu);
500 else if(Enu >= 15 && Enu <= 21) ed15_21GeV.Fill(Efit - Enu);
501
502 } else {
503
504 e0.Fill(volume.getX(Elep, true));
505 er.Fill(volume.getX(Efit) - volume.getX(Elep));
506 ee.Fill(volume.getX(Elep), volume.getX(Efit));
507 ea.Fill(Efit - Elep);
508 hzenith.Fill(Elep, zenith);
509 QE.put(log10(Efit/Elep));
510 herrorE.Fill(x, volume.getX(Efit) - volume.getX(Elep));
511 hfracE.Fill(x, fabs(volume.getX(Efit) - volume.getX(Elep))/volume.getX(Elep));
512
513 if(Elep >= 3 && Elep <= 5) ed3_5GeV.Fill(Efit - Elep);
514 else if(Elep >= 8 && Elep <= 11) ed8_11GeV.Fill(Efit - Elep);
515 else if(Elep >= 15 && Elep <= 21) ed15_21GeV.Fill(Efit - Elep);
516
517 }
518 for (int i = 1; i <= ee.GetNbinsX(); ++i) {
519 std::vector<double> Erecovalues;
520 for (int j = 1; j <= ee.GetNbinsY(); ++j) {
521 double binContent = ee.GetBinContent(i, j);
522 // Add E values based on the bin content (number of events)
523 for (int k = 0; k < binContent; ++k) {
524 double Erecobin = ee.GetYaxis()->GetBinCenter(j);
525 Erecovalues.push_back(Erecobin);
526 }
527 }
528 if(Erecovalues.empty()) continue;
529 eeres.SetBinContent(i, Erecovalues[int(Erecovalues.size()*0.5)]);
530 eeres16.SetBinContent(i, Erecovalues[int(Erecovalues.size()*0.16)]);
531 eeres84.SetBinContent(i, Erecovalues[int(Erecovalues.size()*0.84)]);}
532 e2.Fill(volume.getX(Efit, true));
533
534 h2.Fill(x, beta);
535 }
536 }
537 }
538
539 he.Fill(x, W);}
540
541 STATUS(endl);
542
543 NOTICE("Number of events input " << setw(8) << right << job.GetBinContent(1) << endl);
544 if(!isMuon){
545 NOTICE("Number of events with electron " << setw(8) << right << job.GetBinContent(3) << endl);
546 } else{
547 NOTICE("Number of events with muon " << setw(8) << right << job.GetBinContent(3) << endl);
548 }
549 NOTICE("Number of events with fit " << setw(8) << right << job.GetBinContent(4) << endl);
550 NOTICE("Number of events selected " << setw(8) << right << job.GetBinContent(5) << endl);
551 NOTICE("Number of events with neutrino " << setw(8) << right << job.GetBinContent(6) << endl);
552 NOTICE("Number of events contained " << setw(8) << right << job.GetBinContent(7) << endl);
553
554 if (Q.getCount() != 0) {
555 NOTICE("Median, 70, 80 and 90% quantiles of space angle [deg] " << FIXED (6,3) << Q.getQuantile(0.5) << " "<< Q.getQuantile(0.7) << " "<< Q.getQuantile(0.8)<< " "<< Q.getQuantile(0.9)<< endl);
556 }
557 if (QE.getCount() != 0) {
558 NOTICE("Median, 70%, 80% and 90% quantiles of energy ratio log(Ereco/Etrue) " << FIXED (6,3) << QE.getQuantile(0.5) << " "<< QE.getQuantile(0.7) << " "<< QE.getQuantile(0.8)<< " "<< QE.getQuantile(0.9) << endl);
559 }
560 if(Qp.getCount() != 0){
561 NOTICE("Median distance to vertex: " << setw(8) << Qp.getQuantile(0.5) << " "<< Qp.getQuantile(0.7)<<" "<< Qp.getQuantile(0.8)<< " "<<Qp.getQuantile(0.9)<<endl);
562 NOTICE("Median time to vertex: " << setw(8) << Qt.getQuantile(0.5) << " "<< Qt.getQuantile(0.7)<<" "<< Qt.getQuantile(0.8)<< " "<<Qt.getQuantile(0.9)<<endl);
563 NOTICE("Median 4D-distance to vertex: " << setw(8) << Q4d.getQuantile(0.5) << " "<< Q4d.getQuantile(0.7)<<" "<< Q4d.getQuantile(0.8)<< " "<<Q4d.getQuantile(0.9)<<endl);
564 }
565 out.Write();
566 out.Close();
567}
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
NOTICE
#define NOTICE(A)
Definition JMessage.hh:64
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:69
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 PhysicsEvent.hh:434
JAANET::JHead
Monte Carlo run header.
Definition JHead.hh:1236
JGEOMETRY2D::JCircle2D
Data structure for circle in two dimensions.
Definition JCircle2D.hh:35
JGEOMETRY2D::JPosition2D
Data structure for position in two dimensions.
Definition JPosition2D.hh:33
JGEOMETRY2D::JVector2D
Data structure for vector in two dimensions.
Definition JVector2D.hh:34
JGEOMETRY2D::JVector2D::getLengthSquared
double getLengthSquared() const
Get length squared.
Definition JVector2D.hh:188
JGEOMETRY3D::JAxis3D::getIntersection
double getIntersection(const JVector3D &pos) const
Get longitudinal position along axis of position of closest approach with given position.
Definition JAxis3D.hh:146
JGEOMETRY3D::JCylinder3D::is_inside
bool is_inside(const JVector3D &pos) const
Check whether given point is inside cylinder.
Definition JCylinder3D.hh:208
JGEOMETRY3D::JCylinder3D::add
JCylinder3D & add(const JVector3D &pos)
Add position.
Definition JCylinder3D.hh:156
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition JPosition3D.hh:38
JGEOMETRY3D::JPosition3D::getPosition
const JPosition3D & getPosition() const
Get position.
Definition JPosition3D.hh:130
JGEOMETRY3D::JTrack3D::getT
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition JTrack3D.hh:126
JGEOMETRY3D::JTrack3D::add
JTime & add(const JTime &value)
Addition operator.
Definition JGeometry3D/JTime.hh:67
JGEOMETRY3D::JTrack3D::getDirection
JVersor3D getDirection(const JVector3D &pos) const
Get photon direction of Cherenkov light on PMT.
Definition JTrack3D.hh:147
JGEOMETRY3D::JTrack3D::sub
JTime & sub(const JTime &value)
Subtraction operator.
Definition JGeometry3D/JTime.hh:81
JGEOMETRY3D::JTrack3E
3D track with energy.
Definition JTrack3E.hh:32
JGEOMETRY3D::JVector3D
Data structure for vector in three dimensions.
Definition JVector3D.hh:36
JGEOMETRY3D::JVector3D::add
JVector3D & add(const JVector3D &vector)
Add vector.
Definition JVector3D.hh:142
JGEOMETRY3D::JVector3D::getLength
double getLength() const
Get length.
Definition JVector3D.hh:246
JGEOMETRY3D::JVector3D::getLengthSquared
double getLengthSquared() const
Get length squared.
Definition JVector3D.hh:235
JGEOMETRY3D::JVersor3D::getTheta
double getTheta() const
Get theta angle.
Definition JVersor3D.hh:128
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JPHYSICS::JGeanz::getMaximum
double getMaximum(const double E) const
Get depth of shower maximum.
Definition JGeanz.hh:162
JRECONSTRUCTION::JAtmosphericMuon
Auxiliary class to evaluate atmospheric muon hypothesis.
Definition JReconstruction/JEvtToolkit.hh:735
JRECONSTRUCTION::JPointing
Auxiliary class to compare fit results with respect to a reference direction (e.g....
Definition JReconstruction/JEvtToolkit.hh:550
JRECONSTRUCTION::JPosition
Auxiliary class to compare fit results with respect to a reference position (e.g. true muon).
Definition JReconstruction/JEvtToolkit.hh:646
JTOOLS::JRange
Range of values.
Definition JRange.hh:42
JTOOLS::JRange::getLowerLimit
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
JTOOLS::JRange::getUpperLimit
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
KM3NETDAQ::JDAQEvent::end
const_iterator< T > end() const
Get end of data.
KM3NETDAQ::JDAQEvent::begin
const_iterator< T > begin() const
Get begin of data.
time_converter
Auxiliary class to convert DAQ hit time to/from Monte Carlo hit time.
Definition time_converter.hh:20
JSHOWERPOSITIONFIT
static const int JSHOWERPOSITIONFIT
Definition definitions/reconstruction.hh:24
JSHOWERDIRECTIONPREFIT
static const int JSHOWERDIRECTIONPREFIT
Definition definitions/reconstruction.hh:29
JSHOWERPOINTSIMPLEX
static const int JSHOWERPOINTSIMPLEX
Definition definitions/reconstruction.hh:28
JSHOWERCOMPLETEFIT
static const int JSHOWERCOMPLETEFIT
Definition definitions/reconstruction.hh:25
JSHOWER_BJORKEN_Y
static const int JSHOWER_BJORKEN_Y
Definition definitions/reconstruction.hh:26
JSHOWERPREFIT
static const int JSHOWERPREFIT
Definition definitions/reconstruction.hh:23
JSHOWERENERGYPREFIT
static const int JSHOWERENERGYPREFIT
Definition definitions/reconstruction.hh:27
JAANET::getOrigin
Vec getOrigin(const JHead &header)
Get origin.
Definition JHeadToolkit.hh:309
JAANET::getDirection
JDirection3D getDirection(const Vec &dir)
Get direction.
Definition JAAnetToolkit.hh:214
JAANET::getTrack
JTrack3E getTrack(const Trk &track)
Get track.
Definition JAAnetToolkit.hh:262
JAANET::is_electron
bool is_electron(const Trk &track)
Test whether given track is a (anti-)electron.
Definition JAAnetToolkit.hh:330
JAANET::getCylinder
JCylinder3D getCylinder(const JHead &header)
Get cylinder corresponding to the positions of generated tracks (i.e.
Definition JHeadToolkit.hh:331
JAANET::has_neutrino
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Definition JAAnetToolkit.hh:491
JAANET::getPosition
JPosition3D getPosition(const Vec &pos)
Get position.
Definition JAAnetToolkit.hh:178
JAANET::is_muon
bool is_muon(const Trk &track)
Test whether given track is a (anti-)muon.
Definition JAAnetToolkit.hh:338
JAANET::getOffset
Vec getOffset(const JHead &header)
Get offset.
Definition JHeadToolkit.hh:275
JAANET::get_neutrino
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
Definition JAAnetToolkit.hh:513
JGEOMETRY3D::getAngle
double getAngle(const JQuaternion3D &first, const JQuaternion3D &second)
Get space angle between quanternions.
Definition JQuaternion3D.hh:1165
JLANG::getCount
size_t getCount(const array_type< T > &buffer, const JCompare_t &compare)
Count number of unique values.
Definition JVectorize.hh:261
JMATH::pow
T pow(const T &x, const double y)
Power .
Definition JMath.hh:97
JMATH::PI
static const double PI
Mathematical constants.
Definition JMath/JConstants.hh:20
JPHYSICS::geanz
static const JGeanz geanz(1.85, 0.62, 0.54)
Function object for longitudinal EM-shower profile.
JPHYSICS::getIndexOfRefraction
double getIndexOfRefraction()
Get average index of refraction of water corresponding to group velocity.
Definition JPhysics/JConstants.hh:134
JPHYSICS::getInverseSpeedOfLight
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition JPhysics/JConstants.hh:123
JPHYSICS::getSpeedOfLight
const double getSpeedOfLight()
Get speed of light.
Definition JPhysics/JConstants.hh:112
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:43
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
JSUPPORT::getCommonHeader
Head getCommonHeader(const JMultipleFileScanner_t &file_list)
Get common Monte Carlo header.
Definition JMonteCarloFileSupportkit.hh:559
JTOOLS::j
int j
Definition JPolint.hh:792
KM3NETDAQ
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
create_update_issues.origin
str origin
Definition create_update_issues.py:4
Evt
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition Evt.hh:21
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JAANET::JVolume
Auxiliary class for histogramming of effective volume.
Definition JVolume.hh:29
JACOUSTICS::JEvt
Acoustic event fit.
Definition JAcoustics/JEvt.hh:304
JFIT::JHistory::is_event
Auxiliary class to test history.
Definition JHistory.hh:105
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:288
JSUPPORT::JLimit
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
JSUPPORT::JLimit::getLimit
const JLimit & getLimit() const
Get limit.
Definition JLimit.hh:84
JSUPPORT::JLimit::max
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
JSUPPORT::JTriggeredFileScanner
Auxiliary class to synchronously read DAQ events and Monte Carlo events (and optionally other events)...
Definition JTriggeredFileScanner.hh:41
JTOOLS::JQuantile
Auxiliary data structure for running average, standard deviation and quantiles.
Definition JQuantile.hh:46
Trk
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition Trk.hh:15
Trk::E
double E
Energy [GeV] (either MC truth or reconstructed)
Definition Trk.hh:20
quality_sorter
Reconstruction type dependent comparison of track quality.
Definition tools/reconstruction.hh:85
Generated by doxygen 1.12.0