GnumiTree.cxx

Go to the documentation of this file.

#include "GnumiTree.h"
#include "reweight.h"

#include <TTree.h>
#include <TH2D.h>

#include <iostream>
using namespace std;


void GnumiTree::SetTree(TTree& the_tree)
{
    // Got 'a do it this way because all trees are at top level...
    the_tree.SetBranchAddress("run",&run);
    the_tree.SetBranchAddress("evtno",&evtno);
    the_tree.SetBranchAddress("ndxdz",&ndxdz);
    the_tree.SetBranchAddress("ndydz",&ndydz);
    the_tree.SetBranchAddress("npz",&npz);
    the_tree.SetBranchAddress("nenergy",&nenergy);
    the_tree.SetBranchAddress("ndxdznea",&ndxdznea);
    the_tree.SetBranchAddress("ndydznea",&ndydznea);
    the_tree.SetBranchAddress("nenergyn",&nenergyn);
    the_tree.SetBranchAddress("nwtnear",&nwtnear);
    the_tree.SetBranchAddress("ndxdzfar",&ndxdzfar);
    the_tree.SetBranchAddress("ndydzfar",&ndydzfar);
    the_tree.SetBranchAddress("nenergyf",&nenergyf);
    the_tree.SetBranchAddress("nwtfar",&nwtfar);
    the_tree.SetBranchAddress("norig",&norig);
    the_tree.SetBranchAddress("ndecay",&ndecay);
    the_tree.SetBranchAddress("ntype",&ntype);
    the_tree.SetBranchAddress("vx",&vx);
    the_tree.SetBranchAddress("vy",&vy);
    the_tree.SetBranchAddress("vz",&vz);
    the_tree.SetBranchAddress("pdpx",&pdpx);
    the_tree.SetBranchAddress("pdpy",&pdpy);
    the_tree.SetBranchAddress("pdpz",&pdpz);
    the_tree.SetBranchAddress("ppdxdz",&ppdxdz);
    the_tree.SetBranchAddress("ppdydz",&ppdydz);
    the_tree.SetBranchAddress("pppz",&pppz);
    the_tree.SetBranchAddress("ppenergy",&ppenergy);
    the_tree.SetBranchAddress("ppmedium",&ppmedium);
    the_tree.SetBranchAddress("ptype",&ptype);
    the_tree.SetBranchAddress("ppvx",&ppvx);
    the_tree.SetBranchAddress("ppvy",&ppvy);
    the_tree.SetBranchAddress("ppvz",&ppvz);
    the_tree.SetBranchAddress("muparpx",&muparpx);
    the_tree.SetBranchAddress("muparpy",&muparpy);
    the_tree.SetBranchAddress("muparpz",&muparpz);
    the_tree.SetBranchAddress("mupare",&mupare);
    the_tree.SetBranchAddress("necm",&necm);
    the_tree.SetBranchAddress("nimpwt",&nimpwt);
    the_tree.SetBranchAddress("xpoint",&xpoint);
    the_tree.SetBranchAddress("ypoint",&ypoint);
    the_tree.SetBranchAddress("zpoint",&zpoint);
    the_tree.SetBranchAddress("tvx",&tvx);
    the_tree.SetBranchAddress("tvy",&tvy);
    the_tree.SetBranchAddress("tvz",&tvz);
    the_tree.SetBranchAddress("tpx",&tpx);
    the_tree.SetBranchAddress("tpy",&tpy);
    the_tree.SetBranchAddress("tpz",&tpz);
    the_tree.SetBranchAddress("tptype",&tptype);    

    tree = &the_tree;
}


void GnumiTree::paint_flux(TH2D* hist, double npot, double zloc)
{
        int nbinsx = hist->GetNbinsX();
        int nbinsy = hist->GetNbinsY();

        double xmin = hist->GetXaxis()->GetXmin();
        double xmax = hist->GetXaxis()->GetXmax();

        double ymin = hist->GetYaxis()->GetXmin();
        double ymax = hist->GetYaxis()->GetXmax();

        double xstep = (xmax-xmin)/nbinsx;
        double ystep = (ymax-ymin)/nbinsy;

        int nentries = tree->GetEntries();

        cerr << nentries << " entries " <<
            " x:(" << xmin << " -> " << xmax << " += " << xstep << ")"
            " y:(" << ymin << " -> " << ymax << " += " << ystep << ")"
             << endl;

        for (int ind = 0; ind<nentries; ++ind) {
            tree->GetEntry(ind);

            if (ind%10000 == 1) cerr << " " << ind << ", ";

            int neutrino_type = ntype;

            if (neutrino_type != 56) {
                //cerr << "Got " << neutrino_type << " not 56\n";
                continue;
            }

            double vertex[3] = 
                { vx, vy, vz };
            double neutrino_momentum[3] = 
                { ndxdz*npz, ndydz*npz, npz };
            double parent_momentum[3] = 
                { pdpx, pdpy, pdpz };
            double muparent_momentum[3] = 
                { muparpx, muparpy, muparpz };

            double standard_weight = nimpwt/3.14159;
            int parent_type = ptype;        
            // Paint the different reweights
            for (double x = xmin; x<xmax; x+=xstep) {
                for (double y = ymin; y<ymax; y+=ystep) {
                    double det_pos[] = { x, y, zloc } ;
                    double eneu, weight;
                    reweight_neutrino(det_pos,
                                      vertex,
                                      neutrino_type,
                                      neutrino_momentum,
                                      parent_type,
                                      parent_momentum,
                                      muparent_momentum,
                                      eneu, weight);
                    hist->Fill(x,y,standard_weight*weight/npot);
                }
            }
        }
    }

---