19#include "BDSOutputROOTEventBeam.hh"
20#include "BDSOutputROOTEventHistograms.hh"
21#include "BDSOutputROOTEventLoss.hh"
22#include "BDSOutputROOTEventTrajectory.hh"
25#include "EventAnalysis.hh"
26#include "HistogramMeanFromFile.hh"
27#include "PerEntryHistogramSet.hh"
28#include "RBDSException.hh"
29#include "SamplerAnalysis.hh"
33#include "TDirectory.h"
45 Analysis(
"Event.", nullptr,
"EventHistogramsMerged"),
49 processSamplers(false),
50 emittanceOnTheFly(false),
58 bool perEntryAnalysis,
59 bool processSamplersIn,
62 double printModuloFraction,
63 bool emittanceOnTheFlyIn,
64 long int eventStartIn,
66 const std::string& primaryParticleName):
67 Analysis(
"Event.", chainIn,
"EventHistogramsMerged", perEntryAnalysis, debugIn),
71 processSamplers(processSamplersIn),
72 emittanceOnTheFly(emittanceOnTheFlyIn),
73 eventStart(eventStartIn),
75 nEventsToProcess(eventEndIn - eventStartIn)
105 if (!primaryParticleName.empty())
110 {
throw RBDSException(
"No samplers and no particle name - unable to calculate optics without mass of particle");}
118 std::cout <<
"Analysis on \"" << treeName <<
"\" beginning" << std::endl;
123 TH1::AddDirectory(kTRUE);
124 TH2::AddDirectory(kTRUE);
125 TH3::AddDirectory(kTRUE);
126 BDSBH4DBase::AddDirectory(kTRUE);
128 PreparePerEntryHistogramSets();
133 std::cout <<
"Analysis on \"" << treeName <<
"\" complete" << std::endl;
143EventAnalysis::~EventAnalysis() noexcept
156 {std::cout << __METHOD_NAME__ <<
"Entries: " << chain->GetEntries() <<
" " << std::endl;}
164 <<
") in file(s) -> curtailing to # of entries!" << std::endl;
167 bool firstLoop =
true;
171 {CheckSpectraBranches();}
177 std::cout <<
"\rEvent #" << std::setw(8) << i <<
" of " <<
entries;
181 {std::cout << std::endl;}
192 AccumulatePerEntryHistogramSets(i);
198 std::cout << __METHOD_NAME__ << i << std::endl;
201 std::cout << __METHOD_NAME__ <<
"Vector lengths" << std::endl;
202 std::cout << __METHOD_NAME__ <<
"primaries=" <<
event->Primary->n << std::endl;
203 std::cout << __METHOD_NAME__ <<
"eloss=" <<
event->Eloss->n << std::endl;
204 std::cout << __METHOD_NAME__ <<
"nprimary=" <<
event->PrimaryFirstHit->n << std::endl;
205 std::cout << __METHOD_NAME__ <<
"nlast=" <<
event->PrimaryLastHit->n << std::endl;
206 std::cout << __METHOD_NAME__ <<
"ntunnel=" <<
event->TunnelHit->n << std::endl;
207 std::cout << __METHOD_NAME__ <<
"ntrajectory=" <<
event->Trajectory->n << std::endl;
216 std::cout <<
"\rSampler analysis complete " << std::endl;
219void EventAnalysis::CheckSpectraBranches()
228 TerminatePerEntryHistogramSets();
233 std::vector<double> emittance = {0,0,0,0};
246 auto setDefinitions =
Config::Instance()->EventHistogramSetDefinitionsSimple();
247 for (
auto definition : setDefinitions)
257 std::string peSetsDirName =
"PerEntryHistogramSets";
258 std::string siSetsDirName =
"SimpleHistogramSets";
259 std::string cleanedName = treeName;
260 TDirectory* treeDir = outputFile->GetDirectory(cleanedName.c_str());
261 TDirectory* peSetsDir = treeDir->mkdir(peSetsDirName.c_str());
262 TDirectory* siSetsDir = treeDir->mkdir(siSetsDirName.c_str());
267 {s->Write(peSetsDir);}
274 for (
auto h : set.second)
289 std::vector<double> xOpticsPoint;
290 std::vector<double> yOpticsPoint;
291 std::vector<double> lOpticsPoint;
292 xOpticsPoint.resize(25);
293 yOpticsPoint.resize(25);
294 lOpticsPoint.resize(25);
297 TTree* opticsTree =
new TTree(
"Optics",
"Optics");
298 opticsTree->Branch(
"Emitt_x", &(xOpticsPoint[0]),
"Emitt_x/D");
299 opticsTree->Branch(
"Emitt_y", &(yOpticsPoint[0]),
"Emitt_y/D");
300 opticsTree->Branch(
"Alpha_x", &(xOpticsPoint[1]),
"Alpha_x/D");
301 opticsTree->Branch(
"Alpha_y", &(yOpticsPoint[1]),
"Alpha_y/D");
302 opticsTree->Branch(
"Beta_x", &(xOpticsPoint[2]),
"Beta_x/D");
303 opticsTree->Branch(
"Beta_y", &(yOpticsPoint[2]),
"Beta_y/D");
304 opticsTree->Branch(
"Gamma_x", &(xOpticsPoint[3]),
"Gamma_x/D");
305 opticsTree->Branch(
"Gamma_y", &(yOpticsPoint[3]),
"Gamma_y/D");
306 opticsTree->Branch(
"Disp_x", &(xOpticsPoint[4]),
"Disp_x/D");
307 opticsTree->Branch(
"Disp_y", &(yOpticsPoint[4]),
"Disp_y/D");
308 opticsTree->Branch(
"Disp_xp", &(xOpticsPoint[5]),
"Disp_xp/D");
309 opticsTree->Branch(
"Disp_yp", &(yOpticsPoint[5]),
"Disp_yp/D");
310 opticsTree->Branch(
"Mean_x", &(xOpticsPoint[6]),
"Mean_x/D");
311 opticsTree->Branch(
"Mean_y", &(yOpticsPoint[6]),
"Mean_y/D");
312 opticsTree->Branch(
"Mean_xp", &(xOpticsPoint[7]),
"Mean_xp/D");
313 opticsTree->Branch(
"Mean_yp", &(yOpticsPoint[7]),
"Mean_yp/D");
314 opticsTree->Branch(
"Sigma_x", &(xOpticsPoint[8]),
"Sigma_x/D");
315 opticsTree->Branch(
"Sigma_y", &(yOpticsPoint[8]),
"Sigma_y/D");
316 opticsTree->Branch(
"Sigma_xp",&(xOpticsPoint[9]),
"Sigma_xp/D");
317 opticsTree->Branch(
"Sigma_yp",&(yOpticsPoint[9]),
"Sigma_yp/D");
318 opticsTree->Branch(
"S" ,&(xOpticsPoint[10]),
"S/D");
319 opticsTree->Branch(
"Npart" ,&(xOpticsPoint[11]),
"Npart/D");
321 opticsTree->Branch(
"Sigma_Emitt_x", &(xOpticsPoint[12]),
"Sigma_Emitt_x/D");
322 opticsTree->Branch(
"Sigma_Emitt_y", &(yOpticsPoint[12]),
"Sigma_Emitt_y/D");
323 opticsTree->Branch(
"Sigma_Alpha_x", &(xOpticsPoint[13]),
"Sigma_Alpha_x/D");
324 opticsTree->Branch(
"Sigma_Alpha_y", &(yOpticsPoint[13]),
"Sigma_Alpha_y/D");
325 opticsTree->Branch(
"Sigma_Beta_x", &(xOpticsPoint[14]),
"Sigma_Beta_x/D");
326 opticsTree->Branch(
"Sigma_Beta_y", &(yOpticsPoint[14]),
"Sigma_Beta_y/D");
327 opticsTree->Branch(
"Sigma_Gamma_x", &(xOpticsPoint[15]),
"Sigma_Gamma_x/D");
328 opticsTree->Branch(
"Sigma_Gamma_y", &(yOpticsPoint[15]),
"Sigma_Gamma_y/D");
329 opticsTree->Branch(
"Sigma_Disp_x", &(xOpticsPoint[16]),
"Sigma_Disp_x/D");
330 opticsTree->Branch(
"Sigma_Disp_y", &(yOpticsPoint[16]),
"Sigma_Disp_y/D");
331 opticsTree->Branch(
"Sigma_Disp_xp", &(xOpticsPoint[17]),
"Sigma_Disp_xp/D");
332 opticsTree->Branch(
"Sigma_Disp_yp", &(yOpticsPoint[17]),
"Sigma_Disp_yp/D");
333 opticsTree->Branch(
"Sigma_Mean_x", &(xOpticsPoint[18]),
"Sigma_Mean_x/D");
334 opticsTree->Branch(
"Sigma_Mean_y", &(yOpticsPoint[18]),
"Sigma_Mean_y/D");
335 opticsTree->Branch(
"Sigma_Mean_xp", &(xOpticsPoint[19]),
"Sigma_Mean_xp/D");
336 opticsTree->Branch(
"Sigma_Mean_yp", &(yOpticsPoint[19]),
"Sigma_Mean_yp/D");
337 opticsTree->Branch(
"Sigma_Sigma_x", &(xOpticsPoint[20]),
"Sigma_Sigma_x/D");
338 opticsTree->Branch(
"Sigma_Sigma_y", &(yOpticsPoint[20]),
"Sigma_Sigma_y/D");
339 opticsTree->Branch(
"Sigma_Sigma_xp",&(xOpticsPoint[21]),
"Sigma_Sigma_xp/D");
340 opticsTree->Branch(
"Sigma_Sigma_yp",&(yOpticsPoint[21]),
"Sigma_Sigma_yp/D");
342 opticsTree->Branch(
"Mean_E", &(lOpticsPoint[6]),
"Mean_E/D");
343 opticsTree->Branch(
"Mean_t", &(lOpticsPoint[7]),
"Mean_t/D");
344 opticsTree->Branch(
"Sigma_E", &(lOpticsPoint[8]),
"Sigma_E/D");
345 opticsTree->Branch(
"Sigma_t", &(lOpticsPoint[9]),
"Sigma_t/D");
346 opticsTree->Branch(
"Sigma_Mean_E", &(lOpticsPoint[18]),
"Sigma_Mean_E/D");
347 opticsTree->Branch(
"Sigma_Mean_t", &(lOpticsPoint[19]),
"Sigma_Mean_t/D");
348 opticsTree->Branch(
"Sigma_Sigma_E", &(lOpticsPoint[20]),
"Sigma_Sigma_E/D");
349 opticsTree->Branch(
"Sigma_Sigma_t", &(lOpticsPoint[21]),
"Sigma_Sigma_t/D");
351 opticsTree->Branch(
"xyCorrelationCoefficent", &(xOpticsPoint[24]),
"xyCorrelationCoefficent/D");
355 xOpticsPoint = entry[0];
356 yOpticsPoint = entry[1];
357 lOpticsPoint = entry[2];
368 {s->Process(firstTime);}
381void EventAnalysis::PreparePerEntryHistogramSets()
386 auto setDefinitions = c->EventHistogramSetDefinitionsPerEntry();
387 for (
const auto& def : setDefinitions)
392void EventAnalysis::AccumulatePerEntryHistogramSets(
long int entryNumber)
395 {peSet->AccumulateCurrentEntry(entryNumber);}
398void EventAnalysis::TerminatePerEntryHistogramSets()
401 {peSet->Terminate();}
406 std::vector<TH1*> outputHistograms;
Base class for any TTree analysis.
virtual void SimpleHistograms()
Process histogram definitions from configuration instance.
void PreparePerEntryHistograms()
Create structures necessary for per entry histograms.
long int entries
Number of entries in the chain.
bool debug
Whether debug print out is used or not.
HistogramMeanFromFile * histoSum
Merge of per event stored histograms.
void AccumulatePerEntryHistograms(long int entryNumber)
Accumulate means and variances for per entry histograms.
virtual void Write(TFile *outputFile)
Write rebdsim histograms.
virtual void UserProcess()
Virtual function for user to overload and use. Does nothing by default.
bool perEntry
Whether to analyse each entry in the tree in a for loop or not.
void FillHistogram(HistogramDef *definition, std::vector< TH1 * > *outputHistograms=nullptr)
Create an individual histogram based on a definition.
static Config * Instance(const std::string &fileName="", const std::string &inputFilePath="", const std::string &outputFileName="", const std::string &defaultOutputFileSuffix="_ana")
Singleton accessor.
std::vector< PerEntryHistogramSet * > perEntryHistogramSets
Cache of all per entry histogram sets.
void FillHistogram(HistogramDefSet *definition)
Fill a set of simple histograms across all events.
void Initialise()
Initialise each sampler analysis object in samplerAnalysis.
bool printOut
Whether to print out at all per-event.
Event * event
Event object that data loaded from the file will be loaded into.
bool processSamplers
Whether to process samplers.
long int eventEnd
Event index to end analysis at.
virtual void Process()
Operate on each entry in the event tree.
virtual void SimpleHistograms()
Process histogram definitions from configuration instance.
int printModulo
Cache of print modulo fraction.
std::vector< std::vector< std::vector< double > > > opticalFunctions
Optical functions from all samplers.
void ProcessSamplers(bool firstTime=false)
Process each sampler analysis object.
std::vector< SamplerAnalysis * > samplerAnalyses
Holder for sampler analysis objects.
std::map< HistogramDefSet *, std::vector< TH1 * > > simpleSetHistogramOutputs
Map of simple histograms created per histogram set for writing out.
virtual void Write(TFile *outputFileName)
Write analysis including optical functions to an output file.
long int nEventsToProcess
Difference between start and stop.
virtual void Terminate()
Terminate each individual sampler analysis and append optical functions.
bool emittanceOnTheFly
Whether to calculate emittance fresh at each sampler.
long int eventStart
Event index to start analysis from.
void SetPrintModuloFraction(double fraction)
Set how often to print out information about the event.
BDSOutputROOTEventSampler< double > * GetPrimaries()
Accessor.
BDSOutputROOTEventHistograms * Histos
Local variable ROOT data is mapped to.
bool UsePrimaries() const
Whether there is primary data in the output file.
std::vector< BDSOutputROOTEventSampler< double > * > Samplers
Local variable ROOT data is mapped to.
Specification for a set of histograms.
std::vector< HistogramDef * > definitionsV
Vector version for easy iteration.
Accumulator to merge pre-made per-entry histograms.
void Accumulate(BDSOutputROOTEventHistograms *hNew)
Histogram over a set of integers not number line.
General exception with possible name of object and message.
Analysis routines for an individual sampler.
static void UpdateMass(SamplerAnalysis *s)
Set primary particle mass for optical functions from sampler data.