JPizza.cc File Reference

Example program to verify generator data. More...

#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 "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"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Example program to verify generator data.

Author

mdejong

Definition in file JPizza.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 126 of file JPizza.cc.

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