JPostfit.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile.h"

#include "evt/Head.hh"
#include "evt/Evt.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JDAQ/JDAQEvent.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"


namespace {

  using KM3NETDAQ::JDAQEvent;


  /**
   * Count number of unique identifiers in event.
   *
   * \param  __begin       begin of data
   * \param  __end         end   of data
   * \return               number of unique identifiers
   */
  template<class JHit_t>
  inline int getCount(JDAQEvent::const_iterator<JHit_t> __begin, 
                      JDAQEvent::const_iterator<JHit_t> __end)
  {
    using namespace std;
    using namespace KM3NETDAQ;

    typedef JDAQModuleIdentifier  JType_t;
    //typedef JDAQPMTIdentifier     JType_t;

    set<JType_t> buffer;

    for (JDAQEvent::const_iterator<JHit_t> i = __begin; i != __end; ++i) {
      buffer.insert(JType_t(*i));
    }

    return buffer.size();
  }
}


/**
 * \file
 *
 * Example program to histogram fit results.
 * \author mdejong
 */
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_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;
  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)         = JMUONPREFIT, JMUONSIMPLEX, JMUONGANDALF, JMUONENERGY, JMUONSTART, JMUONPATH; 
    zap['a'] = make_field(atmosphere)          = JAtmosphericMuon(90.0, 90.0);
    zap['O'] = make_field(option)              = "E", "N", "LINE", "LOGE";
    zap['d'] = make_field(debug)               =  2;

    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << 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());

  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 muon = event->mc_trks.end();
    
    for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
        
      if (is_muon(*track)) {

        Emu += track->E;
        
        if (track->E > Emax) {
          muon = track;
          Emax = track->E;
        }
      }
    }

    if (muon == 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(*muon));
    
      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(*muon);
        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));

        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);
  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();
}