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 "TDatabasePDG.h"
#include "TParticlePDG.h"

#include "km3net-dataformat/offline/Evt.hh"
#include "JAAnet/JAAnetToolkit.hh"

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

#include "JTools/JRange.hh"

#include "JPhysics/JGeane.hh"

#include "JSirene/JPythia.hh"

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

namespace {

  using namespace JPP;

  /**
   * Get track energy.\n
   * This method includes ionization and radiative energy losses for muons.
   *
   * \param  trk          track
   * \param  start        starting position of track [m]
   * \return              energy [GeV]
   */
  inline double getE(const Trk& trk, const Vec& start)
  {
    return (is_muon(trk) ? gWater(trk.E, -(trk.pos - start).len()) : trk.E);
  }
  
  
  /**
   * Get track kinetic energy.
   * This method includes ionization and radiative energy losses for muons.
   *
   * \param  trk          track
   * \param  start        starting position of track [m]
   * \return              kinetic energy [GeV]
   */
  inline double getEkin(const Trk& trk, const Vec& start)
  {
    const TDatabasePDG pdb;
    
    const double energy = getE(trk, start);
    const double mass   = pdb.GetParticle(trk.type)->Mass();

    if (energy > mass) {
      
      return sqrt((energy - mass) * (energy + mass));
      
    } else {

      return 0.0;
    }
  }
  
  
  /**
   * 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 (has_neutrino(evt) && 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];

        E1 += getE(track, neutrino.pos);
      }
    }

    return E1;
  }


  /**
   * Get visible energy of final state track.\n
   * This method includes ionization and radiative energy losses for muons.
   *
   * \param track          track
   * \param start          starting position of track [m]
   * \return               visible energy [GeV]
   */
  inline double getVisibleEnergy(const Trk& track, const Vec& start) {

    if (is_muon(track)) {

      const double muMeters = gWater.getX(getE(track, start), MASS_MUON); 

      return muMeters / geanc();
          
    } else {
          
      return pythia(track.type, getEkin(track, start));
    }
  }
  

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

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

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

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

        Evis += getVisibleEnergy(evt.mc_trks[i], neutrino.pos);
      }
    }

    return Evis;
  }


  /**
   * Get visible-energy-weighted mean of the track directions.
   *
   * \param  evt         event
   * \return             visible energy [GeV]
   */
  inline Vec getVisibleEnergyDirection(const Evt& evt)
  {
    Vec EvisDir;
    
    const double EvisTotal = getVisibleEnergy(evt);

    if (has_neutrino(evt) && 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];
        const double EvisFrac = getVisibleEnergy(track, neutrino.pos) / EvisTotal;

        EvisDir += EvisFrac * track.dir;
      }
    }

    return EvisDir.normalize();
  }
  

  /**
   * Get the transverse component of the visible energy of a track\n
   * with respect to a reference direction.
   *
   * \param  track        particle track
   * \param  start        starting position of track [m]
   * \param  refDir       reference direction
   * \return              transverse component of visible energy [GeV]
   */
  inline double getTransverseVisibleEnergy(const Trk& track,
                                           const Vec& start,
                                           const Vec& refDir)
  {
    const double costh = getDirection(refDir).getDot(getDirection(track));
    const double sinth = (1 + costh) * (1 - costh);
    
    return sinth * getVisibleEnergy(track, start);
  }

  
  /**
   * Get the square-root of the summed square transverse visible energy of all final state tracks\n
   * defined with respect to a reference direction.
   *
   * \param  track        particle track
   * \param  refDir       reference direction
   * \return              transverse component of visible energy [GeV]
   */
  inline double getTransverseVisibleEnergy(const Evt& evt, const Vec& refDir)
  {
    double EvisT2 = 0.0;
    
    if (has_neutrino(evt) && 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 double EvisT = getTransverseVisibleEnergy(evt.mc_trks[i], neutrino.dir, refDir);
        
        EvisT2 += EvisT * EvisT;
      }
    }

    return sqrt(EvisT2);
  }

  
  /**
   * Get momentum vector of initial state.\n
   * This method assumes neutrino DIS on a stationary nucleus
   *
   * \param  evt         event
   * \return             energy [GeV]
   */
  inline Vec getP0(const Evt& evt)
  {
    Vec P0;
    
    if (has_neutrino(evt) && is_neutrino(evt.mc_trks[0])) {

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

      P0 = neutrino.E * neutrino.dir;
    }

    return P0;
  }


  /**
   * Get momentum vector of final state.\n
   * This method includes muon energy loss.
   *
   * \param  evt         event
   * \return             final state momentum vector [GeV]
   */
  inline Vec getP1(const Evt& evt)
  {
    Vec P1;
    
    if (has_neutrino(evt) && 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];
        const double trackEk = getEkin(track, neutrino.pos);
        
        P1 += trackEk * track.dir;
      }
    }
    
    return P1;
  }
}


/**
 * \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 ht("ht", "Square-root of summed squared transverse visible energy as fraction of initial state energy",
           120,    0.0,     1.2);

  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);
      const Vec&   EvisDir = getVisibleEnergyDirection (*event);
      const double EvisT   = getTransverseVisibleEnergy(*event, EvisDir);

      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 TDatabasePDG  pdg;
          const TParticlePDG* particle = pdg.GetParticle(track.type);

          cout << setw(32) << left << showpos << particle->GetName() << ' ' << 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/E0);
      ht.Fill(EvisT/E0);
      ha.Fill(getDirection(neutrino).getDot(getDirection(EvisDir)));

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