JPizza.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>

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

#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"

#include "JAAnet/JPDB.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JAAnetToolkit.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"

#include "JTools/JRange.hh"

#include "JGeometry2D/JCircle2D.hh"

#include "JGeometry3D/JCylinder3D.hh"

#include "JSirene/JSireneToolkit.hh"
#include "JSirene/JVisibleEnergyToolkit.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Example program to verify generator data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JMultipleFileScanner<Evt> inputFile;
  JLimit_t&       numberOfEvents = inputFile.getLimit();
  string          outputFile;
  JRange<double>  range;
  int             debug;

  try { 

    JParser<> zap("Example program to verify generator data.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)                                                     = "pizza.root";
    zap['n'] = make_field(numberOfEvents)                                                 = JLimit::max();
    zap['R'] = make_field(range,          "fractional energy and momentum conservation")  = JRange<double>(-0.01, +0.01);
    zap['d'] = make_field(debug)                                                          = 0;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  const JHead& head = getHeader(inputFile);

  const JRange<double>& Erange = head.cut_nu.E;
  
  const JCylinder3D can(JCircle2D(JVector2D(0.0, 0.0), head.can.r),
                        head.can.zmin,
                        head.can.zmax);
  
  TFile out(outputFile.c_str(), "recreate");

  TH1D h0  ("h0", "Energy conservation",
            2000,   -100.0,   100.0);
  
  TH1D h1  ("h1", "Momentum conservation",
            1000,      0.0,   100.0);
  
  TH1D job ("job",  "Particle types",
            10001, -5000.5, +5000.5);
  TH2D hv  ("hv", "Logarithmic visible energy as function of logarithmic initial state energy",
            100, log10(Erange.getLowerLimit()), log10(Erange.getUpperLimit()),    
            100, log10(Erange.getLowerLimit()), log10(Erange.getUpperLimit()));
  TH1D ha  ("ha", "Angle between neutrino-direction and visible-energy-weighted direction",
            200,   -1.0,     1.0);

  TH1D hn  ("hn", "Number of muons per event",
            2001,   -0.5, +2000.5);
  TH1D he  ("he", "Muon energies",
            1000,    0.0,    10.0);
  TH2D hp  ("hp", "Muon positions",
            100,     0.0,   2.0e5,
            200,     0.0,   1.5e3);
  
  while (inputFile.hasNext()) {

    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);

    const Evt* event = inputFile.next();

    for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
      job.Fill((double) track->type);
    }

    if (has_neutrino(*event) && is_neutrino(event->mc_trks[0])) {

      const Trk&   neutrino = event->mc_trks[0];

      const Vec&   P0 = getP0(*event);
      const Vec&   P1 = getP1(*event);

      const double E0 = getE0(*event);
      const double E1 = getE1(*event);

      const double Evis       = getVisibleEnergy      (*event, can);
      const Vec    EvisVector = getVisibleEnergyVector(*event, can);

      if (!range((E0 - E1)/E0) || !range((P0 - P1).len()/P0.len()) || debug >= debug_t) {

        NOTICE(endl << FILL(57,'-') << '\n' << FILL());
        NOTICE("event: " << RIGHT(8) << event->mc_id << RIGHT(67) << "energy [GeV] momentum (x, y, z) [GeV] distance [m]" << endl);

        for (size_t i = 0; i != event->mc_trks.size(); ++i) {
        
          const Trk& track = event->mc_trks[i];

          const string name = (JPDB::getInstance().hasPDG(track.type) ?
                               JPDB::getInstance().getPDG(track.type).name :
                               MAKE_STRING("unknown (type: " << track.type << ")"));

          NOTICE(LEFT(32) << showpos << name << ' ' << FIXED(7,3) << track.E << "    " << FIXED(7,3) << track.E * track.dir << "    " << FIXED(7,3) << (track.pos - neutrino.pos).len() << endl);
        }
        NOTICE(LEFT(32) << "balance:" << ' ' << FIXED(7,3) << E0 - E1 << "    " << FIXED(7,3) << P0 - P1 << endl);
      }

      h0.Fill( E0 - E1 );
      h1.Fill((P0 - P1).len());
            
      hv.Fill(log10(E0), log10(Evis));
      ha.Fill(getDirection(neutrino).getDot(getDirection(EvisVector)));
    }


    int n = 0;

    for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {

      if (is_muon(*track)) {
        ++n;
        he.Fill(log10(track->E));
        hp.Fill(track->pos.x*track->pos.x + track->pos.y*track->pos.y, track->pos.z);
      }
    }

    hn.Fill((Double_t) n);
  }
  
  STATUS(endl); 

  out.Write();
  out.Close();

  return 0;
}