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

  cout.tie(&cerr);

  TFile out(outputFile.c_str(), "recreate");

  TH1D h0 ("h0",  "Fractional energy conservation",
           1001,     -1.0,    +1.0);
  TH1D h1 ("h1",  "Fractional momentum conservation",
           1001,     -1.0,    +1.0);
  TH1D job("job", "Particle types",
           10001, -5000.5, +5000.5);

  TH1D hv("hv", "Visible energy as fraction of initial state energy",
            25,    0.0,     2.5);
  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);


  const JHead& head = getHeader(inputFile);
  
  const JCylinder3D can(JCircle2D(JVector2D(0.0, 0.0), head.can.r),
                        head.can.zmin,
                        head.can.zmax);

  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 Vec Evis = getVisibleEnergy(*event, can);

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

        cout << endl << "--------------------------------------------------------" << endl;
        cout << "event: " << setw(8) << event->mc_id << "                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 JParticle& particle = JPDB::getInstance().getPDG(track.type);

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

      hv.Fill(Evis.len() / E0);
      ha.Fill(getDirection(neutrino).getDot(getDirection(Evis / Evis.len())));

      h0.Fill((E0 - E1)/E0);
      h1.Fill((P0 - P1).len()/P0.len());
    }


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