JShowerPostfit.cc File Reference

#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 "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JTools/JConstants.hh"
#include "JTools/JQuantile.hh"
#include "JTrigger/JTimeConverter.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JFitParameters.hh"
#include "JGizmo/JVolume.hh"
#include "JLang/JPredicate.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 fit results.

Author

mdejong and adomi -> This code is an adaptation of JPostfit to JShowerFit reconstruction

Definition in file JShowerPostfit.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 75 of file JShowerPostfit.cc.

{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  typedef JTriggeredFileScanner<JEvt>                       JTriggeredFileScanner_t;
  typedef JTriggeredFileScanner_t::multi_pointer_type       multi_pointer_type;
 
  JTriggeredFileScanner_t inputFile;
  JLimit_t&               numberOfEvents = inputFile.getLimit();
  string                  outputFile;  
  size_t                  numberOfPrefits;
  JQualitySorter          quality_sorter;
  JAtmosphericMuon        atmosphere;
  double                  Emin_GeV;
  double                  NPE;
  int                     application;
  string                  option;
  bool                    isMuon;
  int                     debug;
 
  try {
 
    JParser<> zap("Example program to histogram fit results.");
 
    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['L'] = make_field(quality_sorter)      = JPARSER::initialised();
    zap['E'] = make_field(Emin_GeV)            = 0.0;
    zap['M'] = make_field(NPE)                 = 0.0;
    zap['A'] = make_field(application)         = JSHOWERPREFIT, JSHOWERPOSITIONFIT, JSHOWERCOMPLETEFIT; 
    zap['a'] = make_field(atmosphere)          = JAtmosphericMuon(90.0, 90.0);
    zap['O'] = make_field(option)              = "E", "N", "LINE", "LOGE";
    zap['I'] = make_field(isMuon);
    zap['d'] = make_field(debug)               =  2;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
  
  cout << "APPLICATION " << application << endl;
 
  JHead head;
 
  try {
    head = getHeader(inputFile);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  const JVolume     volume(head, option != "LINE");
  const JPosition3D center   = get<JPosition3D>(head);
  JCylinder3D       cylinder = get<JCylinder3D>(head);
 
  cylinder.add(center);
 
  NOTICE("Reposition can [m]: " << cylinder << endl);
 
  const double EMIN_GEV = 10e3;  // MUPAGE
 
  TFile out(outputFile.c_str(), "recreate");
 
  TH1D job("job", NULL, 100, 0.5, 100.5);
 
  TH1D hn("hn",  NULL,  250,   -0.5, 249.5);
  TH1D hq("hq",  NULL,  300,    0.0, 600.0);
  TH1D hi("hi",  NULL,  101,   -0.5, 100.5);
  TH1D hv("hv",  NULL,  200,   -6.0,   0.0);
  TH1D h1("h1",  NULL,  200,   -2.0,  +2.0);
  TH1D hc("hc",  NULL,  200,   -1.0,  +1.0);
  TH1D hu("hu",  NULL,  400, -1.0e3, 1.0e3);
 
  TH1D hx("hx",  NULL,  100,   -3.0,  +2.3);    // [log(deg)]
  TH1D hd("hd",  NULL,  100,    0.0,  10.0);    // [m]
  TH1D hz("hz",  NULL,  100, -200.0, 200.0);    // [m]
  TH1D ht("ht",  NULL,  100, -100.0, 100.0);    // [ns]
 
  TH1D e0("e0",  NULL,  100, volume.getXmin(), volume.getXmax());       // [log(E)]
  TH1D e1("e1",  NULL,  100, volume.getXmin(), volume.getXmax());       // [log(E)]
  TH1D e2("e2",  NULL,  100, volume.getXmin(), volume.getXmax());       // [log(E)]
  TH1D er("er",  NULL,  100, -5.0, +5.0);       // [log(E/E)]
  TH2D ee("ee",  NULL,
          40, volume.getXmin(), volume.getXmax(),
          40, volume.getXmin(), volume.getXmax());
 
  TH2D hAngle ("hAngle", ";E_{TRUE};angle [#circ]", 100, volume.getXmin(), volume.getXmax(), 360, 0, 180.0);
 
  const Int_t    ny   = 28800;
  const Double_t ymin =     0.0;                 // [deg]
  const Double_t ymax =   180.0;                 // [deg]
  
  vector<double> X;
 
  if (option.find('E') != string::npos) {
 
    for (double x = volume.getXmin(); x <= volume.getXmax(); x += (volume.getXmax() - volume.getXmin()) / 20) {
      X.push_back(x);
    }
 
  } else {
 
    double x = -0.5;
    
    for ( ; x <=  15.5; x +=  1.0) { X.push_back(x); }
    for ( ; x <=  25.5; x +=  2.0) { X.push_back(x); }
    for ( ; x <=  50.5; x +=  5.0) { X.push_back(x); }
    for ( ; x <= 100.5; x += 10.0) { X.push_back(x); }
    for ( ; x <= 250.5; x += 20.0) { X.push_back(x); }
  }
 
  TH2D     h2("h2", NULL, X.size() - 1, X.data(), ny,    ymin,    ymax);
  TProfile he("he", NULL, X.size() - 1, X.data());
  
  TH2D     ha("ha", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
  TH2D     hb("hb", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
 
  JQuantile Q("Angle", true);
  JQuantile O("Omega", true);
 
 
  while (inputFile.hasNext()) {
 
    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
 
    multi_pointer_type ps = inputFile.next();
 
    JDAQEvent* tev   = ps;
    JEvt*      evt   = ps;
    Evt*       event = ps;
 
    const JTimeConverter converter(*event, *tev);
 
    job.Fill(1.0);
 
    double Emu  = 0.0;
    double Emax = 0.0;
 
    vector<Trk>::const_iterator lepton = event->mc_trks.end();  // can be electron or muon
    //    vector<Trk>::const_iterator electron = event->mc_trks.end();
    
    for (vector<Trk>::const_iterator shower = event->mc_trks.begin(); shower != event->mc_trks.end(); ++shower) {
        
      if (!isMuon && is_electron(*shower)) {
 
        Emu += shower->E;
        
        if (shower->E > Emax) {
          lepton = shower;
          Emax = shower->E;
        }
      } else if(isMuon && is_muon(*shower)){
        
        Emu += shower->E;
        if (shower->E > Emax) {
          lepton = shower;
          Emax = shower->E;
        }
      }
    }
 
    if (lepton == event->mc_trks.end()) {
      continue;
    }
 
    job.Fill(3.0);
 
 
    // abscissa
 
    Double_t x = 0.0;
    
    if (option.find('E') != string::npos)
      x = volume.getX(event->mc_trks[0].E);
    else
      x = getCount(tev->begin<JDAQTriggeredHit>(), tev->end<JDAQTriggeredHit>());
 
      
    // weight for efficiency determination
    
    Double_t W = 0.0;
 
    if (!evt->empty()) {
 
      JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_event(application));
 
      if (evt->begin() == __end) {
        continue;
      }
      
      job.Fill(4.0);
 
      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);
      }
 
      const JPointing pointing(getDirection(*lepton));
    
      JEvt::iterator best  = pointing(evt->begin(), __end);
      const Double_t beta  = pointing.getAngle(*best);
      const double   Efit  = best->getE();
      const double   Eraw  = best->getW(JENERGY_ENERGY, numeric_limits<double>::min());
      const double   mip   = best->getW(JSTART_NPE_MIP, numeric_limits<double>::max());
      //const double   npe   = best->getW(JVETO_NPE,            0.0);
      //const int      count = best->getW(JVETO_NUMBER_OF_HITS, 0.0);
 
      // selection of fit result
      
      bool ok = (Efit >= Emin_GeV  &&
                 mip  >= NPE);
      
      if (ok) {
 
        W = 1.0;
 
        job.Fill(5.0);
 
        hn.Fill((Double_t) best->getNDF());
        hq.Fill(best->getQ());
        hi.Fill((Double_t) distance(evt->begin(), best));
        hc.Fill(best->getDZ());
 
        if (( has_neutrino(*event) && get_neutrino(*event).dir.z >= atmosphere.dot2) ||
            (!has_neutrino(*event) && best->getDZ()              >= atmosphere.dot2)) {
          hu.Fill(atmosphere(evt->begin(), __end));
        }
 
        hx.Fill(max(log10(beta), hx.GetXaxis()->GetXmin()));
 
        Q.put(beta);
        
        JTrack3E ta = getTrack(*lepton);
        JTrack3E tb = getTrack(*best);
        
        ta.add(center);
        tb.sub(converter.putTime());
 
        //ta.move(ta.getIntersection(tb), getSpeedOfLight(), gWater);
        //tb.move(tb.getIntersection(ta), getSpeedOfLight(), gWater);
        static_cast<JTrack3D&>(ta).move(ta.getIntersection(tb), getSpeedOfLight());
        static_cast<JTrack3D&>(tb).move(tb.getIntersection(ta), getSpeedOfLight());
        
        hd.Fill((tb.getPosition() - ta.getPosition()).getLength());
 
        if (has_neutrino(*event)) {
          
          job.Fill(6.0);
 
          JPosition3D vertex = getPosition(get_neutrino(*event));
          
          vertex.add(center);
 
          if (cylinder.is_inside(vertex)) {
 
            job.Fill(7.0);
 
            JTrack3E tc = getTrack(*best);
 
            ha.Fill(best->getX()*best->getX() + best->getY()*best->getY(), best->getZ());
            hz.Fill(tc.getIntersection(vertex));
          }
        }
        
        if (best->getE() >= EMIN_GEV) {
           hb.Fill(best->getX()*best->getX() + best->getY()*best->getY(), best->getZ());
        }
        
        ht.Fill(tb.getT() - ta.getT());
 
        e0.Fill(volume.getX(Emu,  true));
        e1.Fill(volume.getX(Eraw, true));
        e2.Fill(volume.getX(Efit, true));
        er.Fill(volume.getX(Efit) - volume.getX(Emu));
        ee.Fill(volume.getX(Emu), volume.getX(Efit));
        
        hAngle.Fill(Efit, beta);
        h2.Fill(x, beta);
 
        if (best->hasW(JSTART_NPE_MIP)) { 
          h1.Fill(log10(best->getW(JSTART_NPE_MIP)));
        }
 
        if (best->hasW(JVETO_NPE) && best->hasW(JVETO_NUMBER_OF_HITS)) {
 
          const double npe   = best->getW(JVETO_NPE);
          const int    count = best->getW(JVETO_NUMBER_OF_HITS);
          const double pv    = TMath::PoissonI(count, npe);
 
          hv.Fill(max(log10(pv), hv.GetXaxis()->GetXmin()));
        }
      }
    }
 
    he.Fill(x, W);
  }
  STATUS(endl);
 
  NOTICE("Number of events input           " << setw(8) << right << job.GetBinContent(1) << endl);
  if(!isMuon){
    NOTICE("Number of events with electron       " << setw(8) << right << job.GetBinContent(3) << endl);
  } else{
    NOTICE("Number of events with muon       " << setw(8) << right << job.GetBinContent(3) << endl);
  }
  NOTICE("Number of events with fit        " << setw(8) << right << job.GetBinContent(4) << endl);
  NOTICE("Number of events selected        " << setw(8) << right << job.GetBinContent(5) << endl);
  NOTICE("Number of events with neutrino   " << setw(8) << right << job.GetBinContent(6) << endl);
  NOTICE("Number of events contained       " << setw(8) << right << job.GetBinContent(7) << endl);
 
  if (Q.getCount() != 0) {
    NOTICE("Median space angle [deg] " << FIXED     (6,3) << Q.getQuantile(0.5) << endl);
  }
  
  out.Write();
  out.Close();
}