BDSIM
BDSIM is a Geant4 extension toolkit for simulation of particle transport in accelerator beamlines.
Loading...
Searching...
No Matches
EventAnalysis.cc
1/*
2Beam Delivery Simulation (BDSIM) Copyright (C) Royal Holloway,
3University of London 2001 - 2023.
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 printOut(false),
48 printModulo(1),
49 processSamplers(false),
50 emittanceOnTheFly(false),
51 eventStart(0),
52 eventEnd(-1),
53 nEventsToProcess(0)
54{;}
55
57 TChain* chainIn,
58 bool perEntryAnalysis,
59 bool processSamplersIn,
60 bool debugIn,
61 bool printOutIn,
62 double printModuloFraction,
63 bool emittanceOnTheFlyIn,
64 long int eventStartIn,
65 long int eventEndIn,
66 const std::string& primaryParticleName):
67 Analysis("Event.", chainIn, "EventHistogramsMerged", perEntryAnalysis, debugIn),
68 event(eventIn),
69 printOut(printOutIn),
70 printModulo(1),
71 processSamplers(processSamplersIn),
72 emittanceOnTheFly(emittanceOnTheFlyIn),
73 eventStart(eventStartIn),
74 eventEnd(eventEndIn),
75 nEventsToProcess(eventEndIn - eventStartIn)
76{
77 // check we get this right for print out normalisation
78 if (eventEndIn == -1)
79 {nEventsToProcess = (long int)chainIn->GetEntries();}
80
82 {// Create sampler analyses if needed
83 // Analyse the primary sampler in the optics too.
84 SamplerAnalysis* sa = nullptr;
85 SamplerAnalysis* pa = nullptr;
86 if (event->UsePrimaries())
87 {
89 samplerAnalyses.push_back(sa);
90 pa = sa;
91 }
92
93 for (const auto& sampler : event->Samplers)
94 {
95 sa = new SamplerAnalysis(sampler, debug);
96 samplerAnalyses.push_back(sa);
97 }
98 if (!event->UsePrimaries())
99 {
100 if (!samplerAnalyses.empty())
101 {pa = samplerAnalyses[0];}
102 }
103
104 chain->GetEntry(0);
105 if (!primaryParticleName.empty())
106 {SamplerAnalysis::UpdateMass(primaryParticleName);}
107 else if (pa)
109 else
110 {throw RBDSException("No samplers and no particle name - unable to calculate optics without mass of particle");}
111 }
112
113 SetPrintModuloFraction(printModuloFraction);
114}
115
117{
118 std::cout << "Analysis on \"" << treeName << "\" beginning" << std::endl;
120 {
121 // ensure new histograms are added to file
122 // crucial for draw command to work as it identifies the histograms by name
123 TH1::AddDirectory(kTRUE);
124 TH2::AddDirectory(kTRUE);
125 TH3::AddDirectory(kTRUE);
126 BDSBH4DBase::AddDirectory(kTRUE);
128 PreparePerEntryHistogramSets();
129 Process();
130 }
132 Terminate();
133 std::cout << "Analysis on \"" << treeName << "\" complete" << std::endl;
134}
135
137{
138 printModulo = (int)std::ceil((double)nEventsToProcess * fraction);
139 if (printModulo <= 0)
140 {printModulo = 1;}
141}
142
143EventAnalysis::~EventAnalysis() noexcept
144{
145 for (auto& sa : samplerAnalyses)
146 {delete sa;}
147 for (auto& hs : perEntryHistogramSets)
148 {delete hs;}
149}
150
152{
153 Initialise();
154
155 if(debug)
156 {std::cout << __METHOD_NAME__ << "Entries: " << chain->GetEntries() << " " << std::endl;}
157
158 // loop over events
159 if (eventEnd < 0)
160 {eventEnd = entries;}
161 if (eventEnd > entries)
162 {
163 std::cerr << "EventEnd " << eventEnd << " > entries (" << entries
164 << ") in file(s) -> curtailing to # of entries!" << std::endl;
166 }
167 bool firstLoop = true;
168 for (auto i = (Long64_t)eventStart; i < (Long64_t)eventEnd; ++i)
169 {
170 if (firstLoop) // ensure samplers setup for spectra before we load data
171 {CheckSpectraBranches();}
172
173 chain->GetEntry(i);
174 // event analysis feedback
175 if (i % printModulo == 0 && printOut)
176 {
177 std::cout << "\rEvent #" << std::setw(8) << i << " of " << entries;
178 if (!debug)
179 {std::cout.flush();}
180 else
181 {std::cout << std::endl;}
182 }
183
184 // merge histograms stored per event in the output
185 if(firstLoop)
187 else
189
190 // per event histograms
192 AccumulatePerEntryHistogramSets(i);
193
194 UserProcess();
195
196 if(debug)
197 {
198 std::cout << __METHOD_NAME__ << i << std::endl;
199 if (processSamplers)
200 {
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;
208 }
209 }
210
211 if (processSamplers)
212 {ProcessSamplers(firstLoop);}
213 if (firstLoop)
214 {firstLoop = false;} // set to false on first pass of loop
215 }
216 std::cout << "\rSampler analysis complete " << std::endl;
217}
218
219void EventAnalysis::CheckSpectraBranches()
220{
221 for (auto s : perEntryHistogramSets)
222 {s->CheckSampler();}
223}
224
226{
228 TerminatePerEntryHistogramSets();
229
230 if (processSamplers)
231 {
232 //vector of emittance values and errors: emitt_x, emitt_y, err_emitt_x, err_emitt_y
233 std::vector<double> emittance = {0,0,0,0};
234 for (auto& samplerAnalysis : samplerAnalyses)
235 {
236 emittance = samplerAnalysis->Terminate(emittance, !emittanceOnTheFly);
237 opticalFunctions.push_back(samplerAnalysis->GetOpticalFunctions());
238 }
239 }
240}
241
243{
245
246 auto setDefinitions = Config::Instance()->EventHistogramSetDefinitionsSimple();
247 for (auto definition : setDefinitions)
248 {FillHistogram(definition);}
249}
250
251void EventAnalysis::Write(TFile* outputFile)
252{
253 // Write rebdsim histograms:
254 Analysis::Write(outputFile);
255
256 // histogram sets done in this derived class because they only apply to the Event tree
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());
263
264 // per entry histogram sets
265 peSetsDir->cd();
266 for (auto s : perEntryHistogramSets)
267 {s->Write(peSetsDir);}
268 outputFile->cd("/");
269
270 // simple histogram sets
271 siSetsDir->cd();
272 for (const auto& set : simpleSetHistogramOutputs)
273 {
274 for (auto h : set.second)
275 {
276 siSetsDir->Add(h);
277 h->Write();
278 }
279 }
280 outputFile->cd("/");
281
282 // We don't need to write out the optics tree if we didn't process samplers
283 // as there's no possibility of optical data.
284 if (!processSamplers)
285 {return;}
286
287 outputFile->cd("/");
288
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);
295
296 // write optical functions
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");
320
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");
341
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");
350
351 opticsTree->Branch("xyCorrelationCoefficent", &(xOpticsPoint[24]), "xyCorrelationCoefficent/D");
352
353 for(const auto& entry : opticalFunctions)
354 {
355 xOpticsPoint = entry[0];
356 yOpticsPoint = entry[1];
357 lOpticsPoint = entry[2];
358 opticsTree->Fill();
359 }
360 opticsTree->Write();
361}
362
364{
365 if (processSamplers)
366 {
367 for (auto s : samplerAnalyses)
368 {s->Process(firstTime);}
369 }
370}
371
373{
374 if (processSamplers)
375 {
376 for (auto s : samplerAnalyses)
377 {s->Initialise();}
378 }
379}
380
381void EventAnalysis::PreparePerEntryHistogramSets()
382{
383 auto c = Config::Instance();
384 if (c)
385 {
386 auto setDefinitions = c->EventHistogramSetDefinitionsPerEntry();
387 for (const auto& def : setDefinitions)
388 {perEntryHistogramSets.push_back(new PerEntryHistogramSet(def, event, chain));}
389 }
390}
391
392void EventAnalysis::AccumulatePerEntryHistogramSets(long int entryNumber)
393{
394 for (auto& peSet : perEntryHistogramSets)
395 {peSet->AccumulateCurrentEntry(entryNumber);}
396}
397
398void EventAnalysis::TerminatePerEntryHistogramSets()
399{
400 for (auto& peSet : perEntryHistogramSets)
401 {peSet->Terminate();}
402}
403
405{
406 std::vector<TH1*> outputHistograms;
407 for (auto def : definition->definitionsV)
408 {Analysis::FillHistogram(def, &outputHistograms);}
409 simpleSetHistogramOutputs[definition] = outputHistograms;
410}
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:155
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.
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.
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.
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.
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.