Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Plot RAT-PAC MC photon data
/**
* Example script for drawing MCPhoton information from RAT-PAC ROOT files.
*
* Usage:
*
* root[0] .L mcphoton.cpp+
* root[1] plotMCData("myfile.root");
*
* A. Mastbaum <mastbaum@hep.upenn.edu>, June 2015
*/
#include <iostream>
#include <TH1F.h>
#include <TCanvas.h>
#include <RAT/DSReader.hh>
#include <RAT/DS/Root.hh>
#include <RAT/DS/MC.hh>
#include <RAT/DS/MCPMT.hh>
void plotMCData(const char* filename) {
// Canvases and histograms
TCanvas* c1 = new TCanvas();
TH1F* htimes = new TH1F("htimes", ";Time (ns);Counts", 250, 20, 90);
TCanvas* c2 = new TCanvas();
TH1F* hlambda = new TH1F("hlambda", ";Wavelength (nm);Counts", 500, 100, 800);
// Set up RAT ROOT file reader
RAT::DSReader dsReader(filename);
RAT::DS::Root* ds;
// Loop over DS events in the file
while (ds = dsReader.NextEvent()) {
// Get the MC data and loop over MCPMTs
RAT::DS::MC* mc = ds->GetMC();
for (size_t ipmt=0; ipmt<mc->GetMCPMTCount(); ipmt++) {
// Get the MCPMT and loop over MCPhotons
RAT::DS::MCPMT* pmt = mc->GetMCPMT(ipmt);
for (size_t iphoton=0; iphoton<pmt->GetMCPhotonCount(); iphoton++) {
// Get the MCPhoton and extract information
RAT::DS::MCPhoton* photon = pmt->GetMCPhoton(iphoton);
float lambda = photon->GetLambda() * 1e6; // mm -> nm
float time = photon->GetHitTime();
std::cout << lambda << " " << time << std::endl;
// Fill the histograms (for real hits)
if (!photon->IsDarkHit()) {
htimes->Fill(time);
hlambda->Fill(lambda);
}
}
}
}
// Draw
c1->cd();
htimes->Draw();
c1->Update();
c2->cd();
hlambda->Draw();
c2->Update();
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment