Jpp  pmt_effective_area_update
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
JAAPostfit.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 "JTools/JRange.hh"
#include "JPhysics/JConstants.hh"
#include "km3net-dataformat/online/JDAQClock.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JMath/JMathToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JSupport.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 analyse track fit results from Evt formatted data.

Author
bofearraigh

Definition in file JAAPostfit.cc.

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 37 of file JAAPostfit.cc.

38 {
39  using namespace std;
40  using namespace JPP;
41  using namespace KM3NETDAQ;
42 
43  JParallelFileScanner<Evt> inputFile;
44  JLimit_t& numberOfEvents = inputFile.getLimit();
45  string detectorFile;
46  string outputFile;
47  int debug;
48 
49  try {
50 
51  JParser<> zap("Example program to analyse track fit results from Evt formatted data.");
52 
53  zap['f'] = make_field(inputFile);
54  zap['n'] = make_field(numberOfEvents) = JLimit::max();
55  zap['a'] = make_field(detectorFile) = "";
56  zap['o'] = make_field(outputFile) = "histogram.root";
57  zap['d'] = make_field(debug) = 2;
58 
59  zap(argc, argv);
60  }
61  catch(const exception& error) {
62  FATAL(error.what() << endl);
63  }
64 
65  const double rad_to_deg = 180/M_PI;
66  const double R_m = 700.0;
67  const int MAX_OVERLAYS = 200;
68  const double zenith_min = -1.0;
69  const double zenith_max = 1.0;
70 
72 
73  if (detectorFile != "") {
74  try {
75  load(detectorFile, detector);
76  }
77  catch(const JException& error) {
78  FATAL(error);
79  }
80  }
81 
82  const JModuleRouter router(detector);
83 
84  TFile out(outputFile.c_str(), "recreate");
85 
86  TH1D job("job", NULL, 100, 0.5, 100.5);
87  TH1D hz ("hz", "; cos(#theta)" , 21 , zenith_min, zenith_max);
88  TH1D ho ("ho", "; #overlays" , MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
89  TH2D hzo ("hzo", "; cos(#theta) ; #overlays" , 21, zenith_min, zenith_max, MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
90  TH2D hxy ("hxy", "; x position ; y position" , 201, -R_m, R_m, 201, -R_m, R_m);
91  TH1D hq ("hq", "; quality parameter" , 100, -200.0, 800);
92  TH1D hb0 ("hb0", "; log(beta0)" , 60, -2, 3.5);
93  TH1D he ("he", "; log(Ereco [GeV])" , 75, -2, 4);
94  TH2D heo ("heo", "; log(Ereco [GeV]) ; #overlays" , 75, -2, 4, MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
95  TH2D hzq ("hzq", "; cos(#theta); quality" , 21, zenith_min, zenith_max, 50, -10.0, 500.0);
96  TH2D hze ("hze", "; cos(#theta); log(Ereco [GeV])", 21, zenith_min, zenith_max, 75, -2, 4);
97  TH2D hzb0("hzb0", "; cos(#theta); log(beta0)" , 21, zenith_min, zenith_max, 60, -2, 3.5);
98 
99  while (inputFile.hasNext()) {
100 
101  STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
102 
103  const Evt * evt = inputFile.next();
104 
105  if ( has_reconstructed_jppmuon(*evt)) {
106 
108 
109  Trk best_trk = get_best_reconstructed_track<JPP_RECONSTRUCTION_TYPE>(*evt );
110 
112 
113  const double Efit = tb.getE();
114 
115  ho .Fill(evt->overlays);
116  hz .Fill(tb.getDZ());
117  hzo.Fill(tb.getDZ(), evt->overlays) ;
118  hxy.Fill(tb.getX(), tb.getY());
119  hq .Fill(best_trk.lik);
120  hzq.Fill(tb.getDZ(), best_trk.lik );
121 
122  if (best_trk.fitinf[JGANDALF_BETA0_RAD]){
123  hb0.Fill(log10(rad_to_deg* best_trk.fitinf[JGANDALF_BETA0_RAD]) );
124  hzb0.Fill(tb.getDZ(), log10(rad_to_deg* best_trk.fitinf[JGANDALF_BETA0_RAD]) );
125  }
126  he .Fill(log10(Efit));
127  hze.Fill(tb.getDZ(), log10(Efit) );
128  heo.Fill(log10(Efit), evt->overlays );
129  }
130  }
131  STATUS(endl);
132 
133  out.Write();
134  out.Close();
135 }
Utility class to parse command line options.
Definition: JParser.hh:1500
General exception.
Definition: JException.hh:23
JTrack3E getTrack(const Trk &track)
Get track.
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:81
Router for direct addressing of module data in detector data structure.
std::string comment
use as you like
Definition: Trk.hh:33
General purpose class for parallel reading of objects from a single file or multiple files...
string outputFile
bool has_reconstructed_jppmuon(const Evt &evt)
Test whether given event has a track with muon reconstruction.
3D track with energy.
Definition: JTrack3E.hh:30
Auxiliary class for defining the range of iterations of objects.
Definition: JLimit.hh:41
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
std::vector< double > fitinf
place to store additional fit info, for jgandalf, see JFitParameters.hh
Definition: Trk.hh:30
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v2.0.0-15-g59d2e2b https://git.km3net.de/common/km3net-dataformat.
double getY() const
Get y position.
Definition: JVector3D.hh:104
int debug
debug level
Definition: JSirene.cc:63
#define FATAL(A)
Definition: JMessage.hh:67
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
double lik
likelihood or lambda value (for aafit, lambda)
Definition: Trk.hh:22
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
const JLimit & getLimit() const
Get limit.
Definition: JLimit.hh:73
do set_variable DETECTOR_TXT $WORKDIR detector
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:13
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