RAT-PAC setup to look at antiprotons

cylinder.geo

// This file lives in data/puck/
{
name: "GEO",
index: "world",
valid_begin: [0, 0],
valid_end: [0, 0],
mother: "", // world volume has no mother
type: "box",
size: [20000.0, 20000.0, 20000.0], // mm, half-length
material: "cryostat_vacuum",
invisible: 1,
}

{
name: "GEO",
index: "tgt",
valid_begin: [0, 0],
valid_end: [0, 0],
mother: "world",
type: "tube",
r_max: 50.0,
size_z: 100.0,
position: [0.0, 0.0, 0.0],
material: "copper",
color: [0.4, 0.4, 0.6, 0.1],
}

piminus.mac

/rat/db/set DETECTOR experiment "puck"
/rat/db/set DETECTOR geo_file "puck/cylinder.geo"

/run/initialize

/rat/proc count
/rat/procset update 10

/rat/proc user
/rat/procset filename "piminus_angle.root"

/tracking/storeTrajectory 1

/generator/add combo gun:point:uniform
/generator/vtx/set pi- 0 0 -1 80000.0
/generator/pos/set 0 0 500

/run/beamOn 50000

TestUser.cc

// This lives in user/

///////////////////////////////////////////////////////////////////////////////
/// \file TestUser.cc
/// \brief An example User Processor
///
/// This example demonstrates how to write a User Processor. MyUserProc
/// (declared and defined below) just creates a 1-D histogram of the number of
/// photoelectrons in each event, filled it as events are generated and writes
/// it to a file called "numpe.root".
///
/// (Of course, you could easily do this just by writing the events to
/// disk using the outroot processor and making the histogram offline.
/// But then you wouldn't get to see how user processors work!)
///
///////////////////////////////////////////////////////////////////////////////

#include <iostream>
#include <string>
#include <TFile.h>
#include <TH1F.h>
#include <TNtuple.h>
#include <RAT/Processor.hh>
#include <RAT/Log.hh>
#include <RAT/DS/Root.hh>

using namespace RAT;

class MyUserProc : public Processor {
public:
  MyUserProc();
  virtual ~MyUserProc();

  virtual Processor::Result DSEvent(DS::Root *ds);

  virtual void SetS(std::string params, std::string value);

protected:
  TFile *f;
  TNtuple* data;
  std::string fFilename;
};


namespace RAT {
  Processor *construct_user_proc(std::string /*userProcName*/) {
    return new MyUserProc;
  }
} // namespace RAT


// Class definition
MyUserProc::MyUserProc()
    : Processor("user"), f(NULL), data(NULL), fFilename("out.root") {}


MyUserProc::~MyUserProc() {
  if (f) {
    f->cd();
    data->Write();
    f->Close();
    delete f;
  }
}


void MyUserProc::SetS(std::string param, std::string value) {
  if (param == "filename")
    fFilename = value;
  else
    throw ParamUnknown(param);
}


Processor::Result MyUserProc::DSEvent(DS::Root *ds) {
  if (!f) {
    f = TFile::Open(fFilename.c_str(), "RECREATE");
    Log::Assert(f != NULL, "MyUserProc: Unable to open ROOT output file.");
    data = new TNtuple("data", "", "pdg:l:e:p:theta");
  }

  RAT::DS::MC* mc = ds->GetMC();
  for (int i=0; i<mc->GetMCTrackCount(); i++) {
    RAT::DS::MCTrack* track = mc->GetMCTrack(i);
    std::string pname = track->GetParticleName();

    if (pname.find("anti_proton") == std::string::npos) {
      continue;
    }

    TVector3 mz(0, 0, -1);

    for (int j=0; j<track->GetMCTrackStepCount(); j++) {
      RAT::DS::MCTrackStep* step = track->GetMCTrackStep(j);
      if (step->GetVolume() == "tgt" &&
          step->GetProcess() == "Transportation" &&
          //step->GetMomentum().Angle(mz) < 0.279) {
          step->GetEndpoint().Dot(step->GetMomentum()) > 0) {
        // Record particle data
        int pdg = track->GetPDGCode();
        float len = track->GetLength();
        float e = step->GetKE();
        float p = step->GetMomentum().Mag();
        float theta = step->GetMomentum().Angle(mz);
        std::cout << pname << " (" << pdg << "):"
                  << " l=" << len
                  << " e=" << e
                  << " p=" << p
                  << " theta=" << theta
                  << std::endl;
        data->Fill(pdg, len, e, p, theta);
        break;
      }
    }
  }

  return Processor::OK;
}