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"),
48 processSamplers(false),
49 emittanceOnTheFly(false),
56 bool perEntryAnalysis,
57 bool processSamplersIn,
59 double printModuloFraction,
60 bool emittanceOnTheFlyIn,
61 long int eventStartIn,
63 const std::string& primaryParticleName):
64 Analysis(
"Event.", chainIn,
"EventHistogramsMerged", perEntryAnalysis, debugIn),
67 processSamplers(processSamplersIn),
68 emittanceOnTheFly(emittanceOnTheFlyIn),
69 eventStart(eventStartIn),
96 if (!primaryParticleName.empty())
101 {
throw RBDSException(
"No samplers and no particle name - unable to calculate optics without mass of particle");}
109 std::cout <<
"Analysis on \"" << treeName <<
"\" beginning" << std::endl;
114 TH1::AddDirectory(kTRUE);
115 TH2::AddDirectory(kTRUE);
116 TH3::AddDirectory(kTRUE);
117 BDSBH4DBase::AddDirectory(kTRUE);
119 PreparePerEntryHistogramSets();
124 std::cout <<
"Analysis on \"" << treeName <<
"\" complete" << std::endl;
134EventAnalysis::~EventAnalysis() noexcept
147 {std::cout << __METHOD_NAME__ <<
"Entries: " << chain->GetEntries() <<
" " << std::endl;}
155 <<
") in file(s) -> curtailing to # of entries!" << std::endl;
158 bool firstLoop =
true;
162 {CheckSpectraBranches();}
168 std::cout <<
"\rEvent #" << std::setw(8) << i <<
" of " <<
entries;
172 {std::cout << std::endl;}
183 AccumulatePerEntryHistogramSets(i);
189 std::cout << __METHOD_NAME__ << i << std::endl;
192 std::cout << __METHOD_NAME__ <<
"Vector lengths" << std::endl;
193 std::cout << __METHOD_NAME__ <<
"primaries=" <<
event->Primary->n << std::endl;
194 std::cout << __METHOD_NAME__ <<
"eloss=" <<
event->Eloss->n << std::endl;
195 std::cout << __METHOD_NAME__ <<
"nprimary=" <<
event->PrimaryFirstHit->n << std::endl;
196 std::cout << __METHOD_NAME__ <<
"nlast=" <<
event->PrimaryLastHit->n << std::endl;
197 std::cout << __METHOD_NAME__ <<
"ntunnel=" <<
event->TunnelHit->n << std::endl;
198 std::cout << __METHOD_NAME__ <<
"ntrajectory=" <<
event->Trajectory->n << std::endl;
207 std::cout <<
"\rSampler analysis complete " << std::endl;
210void EventAnalysis::CheckSpectraBranches()
219 TerminatePerEntryHistogramSets();
224 std::vector<double> emittance = {0,0,0,0};
237 auto setDefinitions =
Config::Instance()->EventHistogramSetDefinitionsSimple();
238 for (
auto definition : setDefinitions)
248 std::string peSetsDirName =
"PerEntryHistogramSets";
249 std::string siSetsDirName =
"SimpleHistogramSets";
250 std::string cleanedName = treeName;
251 TDirectory* treeDir = outputFile->GetDirectory(cleanedName.c_str());
252 TDirectory* peSetsDir = treeDir->mkdir(peSetsDirName.c_str());
253 TDirectory* siSetsDir = treeDir->mkdir(siSetsDirName.c_str());
258 {s->Write(peSetsDir);}
265 for (
auto h : set.second)
280 std::vector<double> xOpticsPoint;
281 std::vector<double> yOpticsPoint;
282 std::vector<double> lOpticsPoint;
283 xOpticsPoint.resize(25);
284 yOpticsPoint.resize(25);
285 lOpticsPoint.resize(25);
288 TTree* opticsTree =
new TTree(
"Optics",
"Optics");
289 opticsTree->Branch(
"Emitt_x", &(xOpticsPoint[0]),
"Emitt_x/D");
290 opticsTree->Branch(
"Emitt_y", &(yOpticsPoint[0]),
"Emitt_y/D");
291 opticsTree->Branch(
"Alpha_x", &(xOpticsPoint[1]),
"Alpha_x/D");
292 opticsTree->Branch(
"Alpha_y", &(yOpticsPoint[1]),
"Alpha_y/D");
293 opticsTree->Branch(
"Beta_x", &(xOpticsPoint[2]),
"Beta_x/D");
294 opticsTree->Branch(
"Beta_y", &(yOpticsPoint[2]),
"Beta_y/D");
295 opticsTree->Branch(
"Gamma_x", &(xOpticsPoint[3]),
"Gamma_x/D");
296 opticsTree->Branch(
"Gamma_y", &(yOpticsPoint[3]),
"Gamma_y/D");
297 opticsTree->Branch(
"Disp_x", &(xOpticsPoint[4]),
"Disp_x/D");
298 opticsTree->Branch(
"Disp_y", &(yOpticsPoint[4]),
"Disp_y/D");
299 opticsTree->Branch(
"Disp_xp", &(xOpticsPoint[5]),
"Disp_xp/D");
300 opticsTree->Branch(
"Disp_yp", &(yOpticsPoint[5]),
"Disp_yp/D");
301 opticsTree->Branch(
"Mean_x", &(xOpticsPoint[6]),
"Mean_x/D");
302 opticsTree->Branch(
"Mean_y", &(yOpticsPoint[6]),
"Mean_y/D");
303 opticsTree->Branch(
"Mean_xp", &(xOpticsPoint[7]),
"Mean_xp/D");
304 opticsTree->Branch(
"Mean_yp", &(yOpticsPoint[7]),
"Mean_yp/D");
305 opticsTree->Branch(
"Sigma_x", &(xOpticsPoint[8]),
"Sigma_x/D");
306 opticsTree->Branch(
"Sigma_y", &(yOpticsPoint[8]),
"Sigma_y/D");
307 opticsTree->Branch(
"Sigma_xp",&(xOpticsPoint[9]),
"Sigma_xp/D");
308 opticsTree->Branch(
"Sigma_yp",&(yOpticsPoint[9]),
"Sigma_yp/D");
309 opticsTree->Branch(
"S" ,&(xOpticsPoint[10]),
"S/D");
310 opticsTree->Branch(
"Npart" ,&(xOpticsPoint[11]),
"Npart/D");
312 opticsTree->Branch(
"Sigma_Emitt_x", &(xOpticsPoint[12]),
"Sigma_Emitt_x/D");
313 opticsTree->Branch(
"Sigma_Emitt_y", &(yOpticsPoint[12]),
"Sigma_Emitt_y/D");
314 opticsTree->Branch(
"Sigma_Alpha_x", &(xOpticsPoint[13]),
"Sigma_Alpha_x/D");
315 opticsTree->Branch(
"Sigma_Alpha_y", &(yOpticsPoint[13]),
"Sigma_Alpha_y/D");
316 opticsTree->Branch(
"Sigma_Beta_x", &(xOpticsPoint[14]),
"Sigma_Beta_x/D");
317 opticsTree->Branch(
"Sigma_Beta_y", &(yOpticsPoint[14]),
"Sigma_Beta_y/D");
318 opticsTree->Branch(
"Sigma_Gamma_x", &(xOpticsPoint[15]),
"Sigma_Gamma_x/D");
319 opticsTree->Branch(
"Sigma_Gamma_y", &(yOpticsPoint[15]),
"Sigma_Gamma_y/D");
320 opticsTree->Branch(
"Sigma_Disp_x", &(xOpticsPoint[16]),
"Sigma_Disp_x/D");
321 opticsTree->Branch(
"Sigma_Disp_y", &(yOpticsPoint[16]),
"Sigma_Disp_y/D");
322 opticsTree->Branch(
"Sigma_Disp_xp", &(xOpticsPoint[17]),
"Sigma_Disp_xp/D");
323 opticsTree->Branch(
"Sigma_Disp_yp", &(yOpticsPoint[17]),
"Sigma_Disp_yp/D");
324 opticsTree->Branch(
"Sigma_Mean_x", &(xOpticsPoint[18]),
"Sigma_Mean_x/D");
325 opticsTree->Branch(
"Sigma_Mean_y", &(yOpticsPoint[18]),
"Sigma_Mean_y/D");
326 opticsTree->Branch(
"Sigma_Mean_xp", &(xOpticsPoint[19]),
"Sigma_Mean_xp/D");
327 opticsTree->Branch(
"Sigma_Mean_yp", &(yOpticsPoint[19]),
"Sigma_Mean_yp/D");
328 opticsTree->Branch(
"Sigma_Sigma_x", &(xOpticsPoint[20]),
"Sigma_Sigma_x/D");
329 opticsTree->Branch(
"Sigma_Sigma_y", &(yOpticsPoint[20]),
"Sigma_Sigma_y/D");
330 opticsTree->Branch(
"Sigma_Sigma_xp",&(xOpticsPoint[21]),
"Sigma_Sigma_xp/D");
331 opticsTree->Branch(
"Sigma_Sigma_yp",&(yOpticsPoint[21]),
"Sigma_Sigma_yp/D");
333 opticsTree->Branch(
"Mean_E", &(lOpticsPoint[6]),
"Mean_E/D");
334 opticsTree->Branch(
"Mean_t", &(lOpticsPoint[7]),
"Mean_t/D");
335 opticsTree->Branch(
"Sigma_E", &(lOpticsPoint[8]),
"Sigma_E/D");
336 opticsTree->Branch(
"Sigma_t", &(lOpticsPoint[9]),
"Sigma_t/D");
337 opticsTree->Branch(
"Sigma_Mean_E", &(lOpticsPoint[18]),
"Sigma_Mean_E/D");
338 opticsTree->Branch(
"Sigma_Mean_t", &(lOpticsPoint[19]),
"Sigma_Mean_t/D");
339 opticsTree->Branch(
"Sigma_Sigma_E", &(lOpticsPoint[20]),
"Sigma_Sigma_E/D");
340 opticsTree->Branch(
"Sigma_Sigma_t", &(lOpticsPoint[21]),
"Sigma_Sigma_t/D");
342 opticsTree->Branch(
"xyCorrelationCoefficent", &(xOpticsPoint[24]),
"xyCorrelationCoefficent/D");
346 xOpticsPoint = entry[0];
347 yOpticsPoint = entry[1];
348 lOpticsPoint = entry[2];
359 {s->Process(firstTime);}
372void EventAnalysis::PreparePerEntryHistogramSets()
377 auto setDefinitions = c->EventHistogramSetDefinitionsPerEntry();
378 for (
const auto& def : setDefinitions)
383void EventAnalysis::AccumulatePerEntryHistogramSets(
long int entryNumber)
386 {peSet->AccumulateCurrentEntry(entryNumber);}
389void EventAnalysis::TerminatePerEntryHistogramSets()
392 {peSet->Terminate();}
397 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.
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.
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.