BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
EventAnalysis.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2022.
4
5This file is part of BDSIM.
6
7BDSIM is free software: you can redistribute it and/or modify
8it under the terms of the GNU General Public License as published
9by the Free Software Foundation version 3 of the License.
10
11BDSIM is distributed in the hope that it will be useful, but
12WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14GNU General Public License for more details.
15
16You should have received a copy of the GNU General Public License
17along with BDSIM. If not, see <http://www.gnu.org/licenses/>.
18*/
19#include "BDSOutputROOTEventBeam.hh"
20#include "BDSOutputROOTEventHistograms.hh"
21#include "BDSOutputROOTEventLoss.hh"
22#include "BDSOutputROOTEventTrajectory.hh"
23#include "Config.hh"
24#include "Event.hh"
25#include "EventAnalysis.hh"
26#include "HistogramMeanFromFile.hh"
27#include "PerEntryHistogramSet.hh"
28#include "RBDSException.hh"
29#include "SamplerAnalysis.hh"
30#include "rebdsim.hh"
31
32#include "TChain.h"
33#include "TDirectory.h"
34#include "TFile.h"
35
36#include <cmath>
37#include <iomanip>
38#include <iostream>
39#include <string>
40#include <vector>
41
42ClassImp(EventAnalysis)
43
45 Analysis("Event.", nullptr, "EventHistogramsMerged"),
46 event(nullptr),
47 printModulo(1),
48 processSamplers(false),
49 emittanceOnTheFly(false),
50 eventStart(0),
51 eventEnd(-1)
52{;}
53
55 TChain* chainIn,
56 bool perEntryAnalysis,
57 bool processSamplersIn,
58 bool debugIn,
59 double printModuloFraction,
60 bool emittanceOnTheFlyIn,
61 long int eventStartIn,
62 long int eventEndIn,
63 const std::string& primaryParticleName):
64 Analysis("Event.", chainIn, "EventHistogramsMerged", perEntryAnalysis, debugIn),
65 event(eventIn),
66 printModulo(1),
67 processSamplers(processSamplersIn),
68 emittanceOnTheFly(emittanceOnTheFlyIn),
69 eventStart(eventStartIn),
70 eventEnd(eventEndIn)
71{
73 {// Create sampler analyses if needed
74 // Analyse the primary sampler in the optics too.
75 SamplerAnalysis* sa = nullptr;
76 SamplerAnalysis* pa = nullptr;
77 if (event->UsePrimaries())
78 {
80 samplerAnalyses.push_back(sa);
81 pa = sa;
82 }
83
84 for (const auto& sampler : event->Samplers)
85 {
86 sa = new SamplerAnalysis(sampler, debug);
87 samplerAnalyses.push_back(sa);
88 }
89 if (!event->UsePrimaries())
90 {
91 if (!samplerAnalyses.empty())
92 {pa = samplerAnalyses[0];}
93 }
94
95 chain->GetEntry(0);
96 if (!primaryParticleName.empty())
97 {SamplerAnalysis::UpdateMass(primaryParticleName);}
98 else if (pa)
100 else
101 {throw RBDSException("No samplers and no particle name - unable to calculate optics without mass of particle");}
102 }
103
104 SetPrintModuloFraction(printModuloFraction);
105}
106
108{
109 std::cout << "Analysis on \"" << treeName << "\" beginning" << std::endl;
111 {
112 // ensure new histograms are added to file
113 // crucial for draw command to work as it identifies the histograms by name
114 TH1::AddDirectory(kTRUE);
115 TH2::AddDirectory(kTRUE);
116 TH3::AddDirectory(kTRUE);
117 BDSBH4DBase::AddDirectory(kTRUE);
119 PreparePerEntryHistogramSets();
120 Process();
121 }
123 Terminate();
124 std::cout << "Analysis on \"" << treeName << "\" complete" << std::endl;
125}
126
128{
129 printModulo = (int)std::ceil((double)entries * fraction);
130 if (printModulo <= 0)
131 {printModulo = 1;}
132}
133
134EventAnalysis::~EventAnalysis() noexcept
135{
136 for (auto& sa : samplerAnalyses)
137 {delete sa;}
138 for (auto& hs : perEntryHistogramSets)
139 {delete hs;}
140}
141
143{
144 Initialise();
145
146 if(debug)
147 {std::cout << __METHOD_NAME__ << "Entries: " << chain->GetEntries() << " " << std::endl;}
148
149 // loop over events
150 if (eventEnd < 0)
151 {eventEnd = entries;}
152 if (eventEnd > entries)
153 {
154 std::cerr << "EventEnd " << eventEnd << " > entries (" << entries
155 << ") in file(s) -> curtailing to # of entries!" << std::endl;
157 }
158 bool firstLoop = true;
159 for (auto i = (Long64_t)eventStart; i < (Long64_t)eventEnd; ++i)
160 {
161 if (firstLoop) // ensure samplers setup for spectra before we load data
162 {CheckSpectraBranches();}
163
164 chain->GetEntry(i);
165 // event analysis feedback
166 if (i % printModulo == 0)
167 {
168 std::cout << "\rEvent #" << std::setw(8) << i << " of " << entries;
169 if (!debug)
170 {std::cout.flush();}
171 else
172 {std::cout << std::endl;}
173 }
174
175 // merge histograms stored per event in the output
176 if(firstLoop)
178 else
180
181 // per event histograms
183 AccumulatePerEntryHistogramSets(i);
184
185 UserProcess();
186
187 if(debug)
188 {
189 std::cout << __METHOD_NAME__ << i << std::endl;
190 if (processSamplers)
191 {
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;
199 }
200 }
201
202 if (processSamplers)
203 {ProcessSamplers(firstLoop);}
204 if (firstLoop)
205 {firstLoop = false;} // set to false on first pass of loop
206 }
207 std::cout << "\rSampler analysis complete " << std::endl;
208}
209
210void EventAnalysis::CheckSpectraBranches()
211{
212 for (auto s : perEntryHistogramSets)
213 {s->CheckSampler();}
214}
215
217{
219 TerminatePerEntryHistogramSets();
220
221 if (processSamplers)
222 {
223 //vector of emittance values and errors: emitt_x, emitt_y, err_emitt_x, err_emitt_y
224 std::vector<double> emittance = {0,0,0,0};
225 for (auto& samplerAnalysis : samplerAnalyses)
226 {
227 emittance = samplerAnalysis->Terminate(emittance, !emittanceOnTheFly);
228 opticalFunctions.push_back(samplerAnalysis->GetOpticalFunctions());
229 }
230 }
231}
232
234{
236
237 auto setDefinitions = Config::Instance()->EventHistogramSetDefinitionsSimple();
238 for (auto definition : setDefinitions)
239 {FillHistogram(definition);}
240}
241
242void EventAnalysis::Write(TFile* outputFile)
243{
244 // Write rebdsim histograms:
245 Analysis::Write(outputFile);
246
247 // histogram sets done in this derived class because they only apply to the Event tree
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());
254
255 // per entry histogram sets
256 peSetsDir->cd();
257 for (auto s : perEntryHistogramSets)
258 {s->Write(peSetsDir);}
259 outputFile->cd("/");
260
261 // simple histogram sets
262 siSetsDir->cd();
263 for (const auto& set : simpleSetHistogramOutputs)
264 {
265 for (auto h : set.second)
266 {
267 siSetsDir->Add(h);
268 h->Write();
269 }
270 }
271 outputFile->cd("/");
272
273 // We don't need to write out the optics tree if we didn't process samplers
274 // as there's no possibility of optical data.
275 if (!processSamplers)
276 {return;}
277
278 outputFile->cd("/");
279
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);
286
287 // write optical functions
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");
311
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");
332
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");
341
342 opticsTree->Branch("xyCorrelationCoefficent", &(xOpticsPoint[24]), "xyCorrelationCoefficent/D");
343
344 for(const auto& entry : opticalFunctions)
345 {
346 xOpticsPoint = entry[0];
347 yOpticsPoint = entry[1];
348 lOpticsPoint = entry[2];
349 opticsTree->Fill();
350 }
351 opticsTree->Write();
352}
353
355{
356 if (processSamplers)
357 {
358 for (auto s : samplerAnalyses)
359 {s->Process(firstTime);}
360 }
361}
362
364{
365 if (processSamplers)
366 {
367 for (auto s : samplerAnalyses)
368 {s->Initialise();}
369 }
370}
371
372void EventAnalysis::PreparePerEntryHistogramSets()
373{
374 auto c = Config::Instance();
375 if (c)
376 {
377 auto setDefinitions = c->EventHistogramSetDefinitionsPerEntry();
378 for (const auto& def : setDefinitions)
379 {perEntryHistogramSets.push_back(new PerEntryHistogramSet(def, event, chain));}
380 }
381}
382
383void EventAnalysis::AccumulatePerEntryHistogramSets(long int entryNumber)
384{
385 for (auto& peSet : perEntryHistogramSets)
386 {peSet->AccumulateCurrentEntry(entryNumber);}
387}
388
389void EventAnalysis::TerminatePerEntryHistogramSets()
390{
391 for (auto& peSet : perEntryHistogramSets)
392 {peSet->Terminate();}
393}
394
396{
397 std::vector<TH1*> outputHistograms;
398 for (auto def : definition->definitionsV)
399 {Analysis::FillHistogram(def, &outputHistograms);}
400 simpleSetHistogramOutputs[definition] = outputHistograms;
401}
Base class for any TTree analysis.
Definition: Analysis.hh:44
virtual void Terminate()
Definition: Analysis.cc:123
virtual void SimpleHistograms()
Process histogram definitions from configuration instance.
Definition: Analysis.cc:86
void PreparePerEntryHistograms()
Create structures necessary for per entry histograms.
Definition: Analysis.cc:100
long int entries
Number of entries in the chain.
Definition: Analysis.hh:99
bool debug
Whether debug print out is used or not.
Definition: Analysis.hh:98
HistogramMeanFromFile * histoSum
Merge of per event stored histograms.
Definition: Analysis.hh:97
void AccumulatePerEntryHistograms(long int entryNumber)
Accumulate means and variances for per entry histograms.
Definition: Analysis.cc:111
virtual void Write(TFile *outputFile)
Write rebdsim histograms.
Definition: Analysis.cc:131
virtual void UserProcess()
Virtual function for user to overload and use. Does nothing by default.
Definition: Analysis.cc:83
bool perEntry
Whether to analyse each entry in the tree in a for loop or not.
Definition: Analysis.hh:100
void FillHistogram(HistogramDef *definition, std::vector< TH1 * > *outputHistograms=nullptr)
Create an individual histogram based on a definition.
Definition: Analysis.cc:166
static Config * Instance(const std::string &fileName="", const std::string &inputFilePath="", const std::string &outputFileName="", const std::string &defaultOutputFileSuffix="_ana")
Singleton accessor.
Definition: Config.cc:153
Event level analysis.
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.
virtual void Execute()
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.
Event loader.
Definition: Event.hh:50
BDSOutputROOTEventSampler< double > * GetPrimaries()
Accessor.
Definition: Event.hh:61
BDSOutputROOTEventHistograms * Histos
Local variable ROOT data is mapped to.
Definition: Event.hh:146
bool UsePrimaries() const
Whether there is primary data in the output file.
Definition: Event.hh:106
std::vector< BDSOutputROOTEventSampler< double > * > Samplers
Local variable ROOT data is mapped to.
Definition: Event.hh:140
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.