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 "evt/Evt.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JPDB.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JTools/JRange.hh"
#include "JPhysics/JGeane.hh"

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

namespace {

  using namespace JPP;

  /**
   * Particle Data Book.
   */
  const JPDB pdb = JPDB::getInstance();


  /**
   * Get energy of initial state.\n
   * This method includes baryon number conservation.
   *
   * \param  evt         event
   * \return             energy [GeV]
   */
  inline double getE0(const Evt& evt)
  {
    double E0 = 0.0;

    if (!evt.mc_trks.empty() && is_neutrino(evt.mc_trks[0])) {

      const Trk& neutrino = evt.mc_trks[0];
      
      E0 += neutrino.E;

      for (size_t i = 1; i != evt.mc_trks.size(); ++i) {

        const Trk& track = evt.mc_trks[i];

        if (track.type == TRACK_TYPE_NEUTRON    ||
            track.type == TRACK_TYPE_SIGMA_MINUS) {
          E0 += MASS_NEUTRON;
        }

        if (track.type == TRACK_TYPE_PROTON     ||
            track.type == TRACK_TYPE_LAMBDA     ||
            track.type == TRACK_TYPE_SIGMA_PLUS) {
          E0 += MASS_PROTON;
        }

        if (track.type == TRACK_TYPE_ANTINEUTRON) {
          E0 -= MASS_NEUTRON;
        }

        if (track.type == TRACK_TYPE_ANTIPROTON ||
            track.type == TRACK_TYPE_ANTILAMBDA) {
          E0 -= MASS_PROTON;
        }
      }
    }

    return E0;
  }


  /**
   * Get energy of final state.\n
   * This method includes muon energy loss.
   *
   * \param  evt         event
   * \return             energy [GeV]
   */
  inline double getE1(const Evt& evt)
  {
    double E1 = 0.0;

    if (!evt.mc_trks.empty() && is_neutrino(evt.mc_trks[0])) {

      const Trk& neutrino = evt.mc_trks[0];

      for (size_t i = 1; i != evt.mc_trks.size(); ++i) {

        const Trk& track = evt.mc_trks[i];

        if (track.type == TRACK_TYPE_MUON       ||
            track.type == TRACK_TYPE_ANTIMUON) {

          E1 += gWater(track.E, -(track.pos - neutrino.pos).len());

        } else {

          E1 += track.E;
        }
      }
    }

    return E1;
  }
}


/**
 * \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 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",  NULL,  1001,    -1.0,    +1.0);
  TH1D job("job", NULL, 10001, -5000.5, +5000.5);

  TH1D hn ("hn",  NULL,  2001,    -0.5, +2000.5);
  TH1D he ("he",  NULL,  1000,     0.0,    10.0); 

  TH2D h2 ("h2",  NULL,
           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 double E0 = getE0(*event);
      const double E1 = getE1(*event);

      if (!range((E0 - E1)/E0) || debug >= debug_t) {

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

        cout << endl << "--------------------------------------------------------" << endl;
        cout << "event: " << setw(8) << event->mc_id << "                energy [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 = pdb.getPDG(track.type);

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

      h0.Fill((E0 - E1)/E0);
    }


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