Jpp - the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JAAPostfit.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 #include "JTools/JRange.hh"
16 #include "JPhysics/JConstants.hh"
18 #include "JAAnet/JAAnetToolkit.hh"
19 #include "JMath/JMathToolkit.hh"
20 #include "JDetector/JDetector.hh"
22 
25 #include "JSupport/JSupport.hh"
26 
27 #include "Jeep/JParser.hh"
28 #include "Jeep/JMessage.hh"
29 
30 
31 /**
32  * \file
33  *
34  * Example program to analyse track fit results from AAnet formatted data.
35  * \author bofearraigh
36  */
37 int main(int argc, char **argv)
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 AAnet 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
int main(int argc, char *argv[])
Definition: Main.cc:15
ROOT TTree parameter settings of various packages.
Auxiliary methods for geometrical methods.
JTrack3E getTrack(const Trk &track)
Get track.
Synchronously read DAQ events and Monte Carlo events (and optionally other events).
#define STATUS(A)
Definition: JMessage.hh:63
Detector data structure.
Definition: JDetector.hh:80
Auxiliary methods for selection of reconstructed tracks.
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
Data structure for detector geometry and calibration.
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
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:30
static const int JGANDALF_BETA0_RAD
KM3NeT Data Definitions v2.0.0-15-g59d2e2b https://git.km3net.de/common/km3net-dataformat.
Physics constants.
double getY() const
Get y position.
Definition: JVector3D.hh:104
int debug
debug level
Definition: JSirene.cc:63
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...
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
Auxiliary class to define a range between two values.
double lik
likelihood or lambda value (for aafit, lambda)
Definition: Trk.hh:22
Utility class to parse command line options.
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