Skip to content

Instantly share code, notes, and snippets.

@mastbaum
Created Apr 11, 2017
Embed
What would you like to do?
RAT-PAC setup to look at antiprotons
// 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],
}
/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
// 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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment