NeutronFates.cxx

Go to the documentation of this file.

// Simple record of the simulated fate of a neutron
//
// Chris Backhouse - c.backhouse@ucl.ac.uk

#include "g4nova/NeutronFates.h"

#include "TFile.h"
#include "TTree.h"

#include <algorithm>

namespace g4n
{
  static_assert(sizeof(NeutronFate::XYZ) == 12, "sizes");
  static_assert(sizeof(NeutronFate::Child) == 32, "sizes");
  static_assert(sizeof(NeutronFate) == 16, "sizes");

  void LoadNeutronFates(const std::string& name,
                        std::vector<NeutronFate>& out,
                        NeutronFate::Child*& childArena)
  {
    TFile fin(name.c_str());
    assert(!fin.IsZombie());
    TTree* tr = (TTree*)fin.Get("tr");
    assert(tr);

    out.reserve(tr->GetEntries());

    NeutronFate nf;
    tr->SetBranchAddress("nchildren", &nf.nchildren);

    long totchildren = 0;
    for(long i = 0; i < tr->GetEntries(); ++i){
      tr->GetEntry(i);
      totchildren += nf.nchildren;
    }
    childArena = new NeutronFate::Child[totchildren];
    NeutronFate::Child* arenaPtr = childArena;

    tr->SetBranchAddress("p", &nf.p);
    //    tr->SetBranchAddress("dt", &nf.dt);
    //    tr->SetBranchAddress("stepz", &nf.stepz);
    //    tr->SetBranchAddress("stepx", &nf.stepx);

    int childpdg[10000];
    tr->SetBranchAddress("childpdg", &childpdg);
    float childx[10000], childy[10000], childz[10000], childt[10000];
    tr->SetBranchAddress("childx", &childx);
    tr->SetBranchAddress("childy", &childy);
    tr->SetBranchAddress("childz", &childz);
    tr->SetBranchAddress("childt", &childt);
    float childpx[10000], childpy[10000], childpz[10000];
    tr->SetBranchAddress("childpx", &childpx);
    tr->SetBranchAddress("childpy", &childpy);
    tr->SetBranchAddress("childpz", &childpz);

    for(long i = 0; i < tr->GetEntries(); ++i){
      tr->GetEntry(i);

      if(nf.nchildren == 0){
        nf.children = 0;
      }
      else{
        nf.children = arenaPtr;
        for(int j = 0; j < nf.nchildren; ++j){
          arenaPtr->dr = {childx[j], childy[j], childz[j]};
          arenaPtr->dt = childt[j];
          arenaPtr->p = {childpx[j], childpy[j], childpz[j]};
          arenaPtr->pdg = childpdg[j];
          ++arenaPtr;
        }
      }

      out.push_back(nf);
    }

    assert(arenaPtr == childArena+totchildren);

    std::sort(out.begin(), out.end(), [](auto a, auto b){return a.p < b.p;});
  }
}