#include <string>
#include <iostream>
#include <iomanip>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.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 "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JVolume.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
/**
* \file
*
* Program to histogram fit results from reconstructed data
* \author bofearraigh, rgracia
*/
int main(int argc, char **argv)
{
using namespace std;
using namespace JPP;
using namespace KM3NETDAQ;
typedef JTriggeredFileScanner<JEvt> JTriggeredFileScanner_t;
typedef JTriggeredFileScanner_t::multi_pointer_type multi_pointer_type;
JTriggeredFileScanner<JEvt> inputFile;
JLimit_t& numberOfEvents = inputFile.getLimit();
string outputFile;
string option;
size_t numberOfPrefits; // The number of reconstructed tracks to plot per event (1 by default)
JQualitySorter quality_sorter; // The criterion by which to pick the "best" reconstructed track for each event.
int application;
JRange<double> energy_range;
JRange<double> coszenith_range;
int debug;
try {
JParser<> zap("Program to histogram fit results from reconstructed data.");
zap['f'] = make_field(inputFile);
zap['o'] = make_field(outputFile) = "postfit.root";
zap['n'] = make_field(numberOfEvents) = JLimit::max();
zap['N'] = make_field(numberOfPrefits) = 1;
zap['A'] = make_field(application) = JMUONENERGY, JMUONPREFIT, JMUONSIMPLEX, JMUONGANDALF, JMUONSTART; //makes histograms after this stage of the reconstruction
zap['L'] = make_field(quality_sorter) = JPARSER::initialised();
zap['E'] = make_field(energy_range) = JRange<double>();
zap['c'] = make_field(coszenith_range) = JRange<double>(-1.0,1.0);
zap['d'] = make_field(debug) = 2;
zap(argc, argv);
}
catch(const JException& error) {
FATAL(error);
}
const double rad_to_deg = 180/M_PI;
const double R_m = 700.0;
const int MAX_OVERLAYS = 300;
const double MAX_TRIGGERED_HITS = 2500;
const double E_RANGE_MIN = -1.0;
const double E_RANGE_MAX = 7;
const double Z_RANGE_MIN = coszenith_range.getLowerLimit();
const double Z_RANGE_MAX = coszenith_range.getUpperLimit();
TFile out(outputFile.c_str(), "recreate");
TH1D job("job", NULL, 100, 0.5, 100.5);
TH1D hz ("hz", "; cos(#theta_{zenith})" , 21 , Z_RANGE_MIN, Z_RANGE_MAX);
TH1D ho ("ho", "; # of overlays" , MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
TH2D hzo ("hzo", "; cos(#theta_{zenith}) ; # of overlays" , 21, Z_RANGE_MIN, Z_RANGE_MAX, MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
TH2D hxy ("hxy", "; x position ; y position" , 201, -R_m, R_m, 201, -R_m, R_m);
TH1D hn ("hn", "; # of triggered hits" , 600, -0.5, MAX_TRIGGERED_HITS-0.5);
TH1D hN ("hN", "; JEnergy # of hits" , 600, -0.5, MAX_TRIGGERED_HITS-0.5);
TH1D hq ("hq", "; quality parameter" , 200, -200.0, 800);
TH1D hb0 ("hb0", "; log(beta0)" , 60, -1, 3.5);
TH1D he ("he", "; log(Ereco [GeV])" , 75, E_RANGE_MIN, E_RANGE_MAX);
TH2D heo ("heo", "; log(Ereco [GeV]) ; # of overlays" , 75, E_RANGE_MIN, E_RANGE_MAX, MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
TH2D hen ("hen", "; log(Ereco [GeV]) ; # of triggered hits" , 75, E_RANGE_MIN, E_RANGE_MAX, 600 , -0.5, MAX_TRIGGERED_HITS-0.5);
TH2D heN ("heN", "; log(Ereco [GeV]) ; JEnergy # of hits" , 75, E_RANGE_MIN, E_RANGE_MAX, 600 , -0.5, MAX_TRIGGERED_HITS-0.5);
TH2D hzq ("hzq", "; cos(#theta_{zenith}); quality" , 21, Z_RANGE_MIN, Z_RANGE_MAX, 600, -200.0, 800.0);
TH2D hzn ("hzn", "; cos(#theta_{zenith}); # of triggered hits", 21, Z_RANGE_MIN, Z_RANGE_MAX, 600 , -0.5, MAX_TRIGGERED_HITS-0.5);
TH2D hzN ("hzN", "; cos(#theta_{zenith}); JEnergy # of hits" , 21, Z_RANGE_MIN, Z_RANGE_MAX, 600, -0.5, MAX_TRIGGERED_HITS-0.5);
TH2D hze ("hze", "; cos(#theta_{zenith}); log(Ereco [GeV])" , 21, Z_RANGE_MIN, Z_RANGE_MAX, 75, E_RANGE_MIN, E_RANGE_MAX);
TH2D hzb0("hzb0", "; cos(#theta_{zenith}); log(beta0)" , 21, Z_RANGE_MIN, Z_RANGE_MAX, 60, E_RANGE_MIN, 3.5);
while (inputFile.hasNext()) { // Loop over each reconstructed event
STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); // if debug > set value, print out this line
multi_pointer_type ps = inputFile.next();
JDAQEvent* tev = ps;
JEvt* evt = ps;
job.Fill(1.0);
if (!evt->empty()) {
JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_event(application));
if (evt->begin() == __end) {
continue;
}
if (numberOfPrefits > 0) {
JEvt::iterator __q = __end;
advance(__end = evt->begin(), min(numberOfPrefits, (size_t) distance(evt->begin(), __q)));
partial_sort(evt->begin(), __end, __q, quality_sorter);
} else {
sort(evt->begin(), __end, quality_sorter);
}
if (numberOfPrefits > 0) {
advance(__end = evt->begin(), min(numberOfPrefits, evt->size()));
}
for (JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) {
JTrack3D track = getTrack(*fit);
if (track.getDZ() > Z_RANGE_MIN && track.getDZ() < Z_RANGE_MAX && fit->getE() > energy_range.getLowerLimit() && fit->getE() < energy_range.getUpperLimit() ) {
const double overlays = tev->getOverlays();
const double Efit = fit->getE();
ho .Fill(overlays);
hn .Fill((double) tev->size<JDAQTriggeredHit>());
hz .Fill(track.getDZ());
hzn.Fill(track.getDZ(), (double) tev->size<JDAQTriggeredHit>());
hzo.Fill(track.getDZ(), overlays) ;
hxy.Fill(track.getX(), track.getY());
hq .Fill(fit->getQ());
hzq.Fill(track.getDZ(), fit->getQ() );
if (fit->hasW(JENERGY_NUMBER_OF_HITS)) {
hN .Fill( fit->getW(JENERGY_NUMBER_OF_HITS) );
heN .Fill(log10(Efit), fit->getW(JENERGY_NUMBER_OF_HITS) );
hzN .Fill(track.getDZ(), fit->getW(JENERGY_NUMBER_OF_HITS) );
}
if (fit->hasW(JGANDALF_BETA0_RAD)) {
hb0 .Fill(log10( rad_to_deg* fit->getW(JGANDALF_BETA0_RAD) ));
hzb0 .Fill(track.getDZ(), log10( rad_to_deg* fit->getW(JGANDALF_BETA0_RAD)) );
}
he .Fill(log10(Efit));
hen .Fill(log10(Efit), (double) tev->size<JDAQTriggeredHit>());
hze .Fill(track.getDZ(), log10(Efit) );
heo .Fill(log10(Efit), overlays );
}
}
}
}
STATUS("Number of events input " << setw(8) << right << job.GetBinContent(1) << endl);
out.Write();
out.Close();
}