19#include "BDSBunchUserFile.hh"
21#include "BDSException.hh"
22#include "BDSGlobalConstants.hh"
23#include "BDSIonDefinition.hh"
24#include "BDSParticleDefinition.hh"
25#include "BDSUtilities.hh"
26#include "BDSWarning.hh"
28#include "parser/beam.h"
30#include "G4IonTable.hh"
31#include "G4ParticleDefinition.hh"
32#include "G4ParticleTable.hh"
36#include "src-external/gzstream/gzstream.h"
55 printedOutFirstTime(false),
56 anEnergyCoordinateInUse(false),
57 changingParticleType(false),
58 endOfFileReached(false),
59 matchDistrFileLength(false)
62 comment = std::regex(
"^\\s*\\#|\\!.*");
69 if (distrFile.empty())
70 {
throw BDSException(
"BDSBunchUserFile::CheckParameters",
"No input file specified for userfile distribution");}
82 endOfFileReached =
false;
83 if (!printedOutFirstTime)
85 G4cout <<
"BDSBunchUserFile::OpenBunchFile> opening " << distrFilePath << G4endl;
86 printedOutFirstTime =
true;
89 InputBunchFile.open(distrFilePath);
90 if (!InputBunchFile.good())
91 {
throw BDSException(
"BDSBunchUserFile::OpenBunchFile>",
"Cannot open bunch file " + distrFilePath);}
97 InputBunchFile.clear();
98 InputBunchFile.close();
99 endOfFileReached =
true;
107 std::regex colons(
":");
108 std::vector<std::string> results;
110 std::sregex_token_iterator end;
111 for (; iter != end; ++iter)
113 std::string res = (*iter).str();
114 results.push_back(res);
116 for (
auto const& token : results)
120 if(token.substr(0,1)==
"E" || token.substr(0,1)==
"P")
123 if (token.substr(0,2)==
"Ek")
125 vari = token.substr(0,2);
126 rest = token.substr(2);
130 vari = token.substr(0,1);
131 rest = token.substr(1);
135 else if (token.substr(0,1)==
"t")
137 vari = token.substr(0, 1);
138 rest = token.substr(1);
141 else if ( ( token.substr(0,1)==
"x" && token.substr(1,1)!=
"p" ) ||
142 ( token.substr(0,1)==
"y" && token.substr(1,1)!=
"p" ) ||
143 ( token.substr(0,1)==
"z" && token.substr(1,1)!=
"p" ) )
145 vari = token.substr(0,1);
146 rest = token.substr(1);
149 else if ( (token.substr(0,2)==
"xp") ||
150 (token.substr(0,2)==
"yp") ||
151 (token.substr(0,2)==
"zp") )
153 vari = token.substr(0,2);
154 rest = token.substr(2);
157 else if (token.substr(0, 1) ==
"S")
159 vari = token.substr(0, 1);
160 rest = token.substr(1);
164 else if(token.substr(0,5)==
"pdgid")
171 else if(token.substr(0,1)==
"w")
177 else if (token.substr(0,1)==
"-")
186 G4String message =
"Cannot determine bunch data format. Failed at token: " + token;
191 std::set<G4String> energyKeys = {
"E",
"Ek",
"P"};
194 std::set<G4String> zKeys = {
"z",
"S"};
201 auto inSet = [=](
const G4String& v){
return s.find(v) != s.end();};
203 int count = std::count_if(
fields.begin(),
fields.end(), inSetD);
206 G4String message =
"More than one of the following set in user file columns (\"distrFileFormat\") ";
207 for (
const auto& st : s)
212 message +=
"\nPossibly conflicting information. Ensure only one by skipping others with \"-\" symbol";
213 throw BDSException(
"BDSBunchUserFile::CheckConflictingParameters>", message);
220 const G4String& rest,
232 G4int pos1 = (G4int)rest.find(
"[");
233 G4int pos2 = (G4int)rest.find(
"]");
234 if (pos1 < 0 || pos2 < 0)
236 G4String message =
"Missing bracket [] in units of \"distrFileFormat\"\n";
237 message +=
"variable : \"" + uName +
"\" and unit \"" + rest +
"\"";
238 throw BDSException(
"BDSBunchUserFile::CheckAndParseUnits>", message);
242 G4String fmt = rest.substr(pos1 + 1, pos2 - 1);
244 sd.unit = unitParser(fmt);
253 for (G4int i = 0; i < nValues; i++)
263 {G4cout <<
"BDSBunchUserFile> ignoring " <<
nlinesIgnore <<
" lines" << G4endl;}
272 G4String msg =
"end of file reached after line " + std::to_string(
lineCounter);
273 msg +=
" before nlinesIgnore ("+std::to_string(
nlinesIgnore)+
") was reached.";
274 throw BDSException(
"BDSBunchUserFile::SkipNLinesIgnoreIntoFile>", msg);
286 {G4cout <<
"BDSBunchUserFile> skipping " <<
nlinesSkip <<
" valid lines" << G4endl;}
292 G4int nLinesValidRead = 0;
300 IncrementNEventsInFileSkipped((
unsigned long long int)
nlinesSkip);
308 G4Transform3D beamlineTransformIn,
309 const G4double beamlineSIn)
325 return std::all_of(line.begin(), line.end(), isspace) || std::regex_search(line, comment);
335 G4long nLinesValid = 0;
353 G4String msg =
"nlinesSkip is greater than the number of valid lines (" + std::to_string(nLinesValidData);
355 throw BDSException(
"BDSBunchUserFile::Initialise>", msg);
359 G4bool nGenerateHasBeenSet = g->NGenerateSet();
360 G4int nEventsPerLoop = (G4int)(nLinesValidData -
nlinesSkip);
362 G4int nAvailable = nEventsPerLoop * distrFileLoopNTimes;
363 G4int nGenerate = g->NGenerate();
364 if (matchDistrFileLength)
366 if (!nGenerateHasBeenSet)
368 g->SetNumberToGenerate(nAvailable);
369 G4cout <<
"BDSBunchUserFile::Initialise> distrFileMatchLength is true -> simulating " << nEventsPerLoop <<
" events";
370 if (distrFileLoopNTimes > 1)
371 {G4cout <<
" " << distrFileLoopNTimes <<
" times";}
376 G4int nEventsRemaining = nAvailable - g->StartFromEvent();
377 g->SetNumberToGenerate(nEventsRemaining);
378 G4cout <<
"BDSBunchUserFile::Initialise> distrFileMatchLength + recreation -> simulate the "
379 << nEventsRemaining <<
" lines left given startFromEvent including possible looping" << G4endl;
384 G4cout <<
"BDSBunchUserFile::Initialise> matchDistrFileLength has been requested "
385 <<
"but ngenerate has been specified -> use ngenerate" << G4endl;
387 if (nGenerate > nAvailable)
389 G4String msg =
"ngenerate (" + std::to_string(nGenerate) +
") is greater than the number of valid lines (";
390 msg += std::to_string(nLinesValidData) +
") and distrFileMatchLength is on.\nChange ngenerate to <= # lines";
391 msg +=
", or don't specifcy ngenerate.\n";
392 msg +=
"This includes nlinesSkip.";
393 throw BDSException(
"BDSBunchUserFile::Initialise>", msg);
399 if ((nGenerate > nEventsPerLoop) && !distrFileLoop)
401 G4String msg =
"ngenerate (" + std::to_string(nGenerate) +
") is greater than the number of valid lines (";
402 msg += std::to_string(nLinesValidData) +
") but distrFileLoop is false in the beam command";
403 throw BDSException(
"BDSBunchUserFile::Initialise>", msg);
416 G4cout <<
"BDSBunchUserFile> End of file reached." << G4endl;
420 G4cout <<
"BDSBunchUserFile> Returning to beginning of file (including nlinesIgnore & nlinesSkip) for next event." << G4endl;
426 {
throw BDSException(__METHOD_NAME__,
"distrFileLoop off but requesting another set of coordinates.");}
433 G4cout <<
"BDSBunchUserFile::RecreateAdvanceToEvent> Advancing file to event: " << eventOffset << G4endl;
434 G4int nEventsPerLoop = nLinesValidData -
nlinesSkip;
435 G4int nAvailable = nEventsPerLoop * distrFileLoopNTimes;
436 G4int nEventsRemaining = nAvailable - eventOffset;
437 if (eventOffset > nEventsPerLoop)
443 G4int nToRoll = eventOffset % nEventsPerLoop;
444 eventOffset = nToRoll;
448 G4String msg =
"eventOffset (" + std::to_string(eventOffset) +
") is greater than the number of valid data lines in this file.\n";
449 msg +=
"This includes nlinesSkip.";
450 throw BDSException(
"BDSBunchUserFile::RecreateAdvanceToEvent>", msg);
455 if (nGenerate > nEventsRemaining && !distrFileLoop)
457 G4String msg =
"ngenerate (" + std::to_string(nGenerate) +
") requested in recreate mode is greater than number\n";
458 msg +=
"of remaining valid lines in file (" + std::to_string(nEventsRemaining) +
") and distrFileLoop is turned off.";
468 for (G4int i = 0; i < eventOffset; i++)
490 G4double E = 0, Ek = 0, P = 0, x = 0, y = 0, z = 0, xp = 0, yp = 0, zp = 0, t = 0;
494 G4bool zpdef =
false;
498 G4bool updateParticleDefinition =
false;
506 G4bool lineIsBad =
true;
527 std::vector<std::string> results;
528 std::regex wspace(
"\\s+");
530 std::sregex_token_iterator iter(line.begin(), line.end(), wspace, -1);
531 std::sregex_token_iterator end;
532 for (; iter != end; ++iter)
534 std::string res = (*iter).str();
535 results.push_back(res);
538 if (results.size() <
fields.size())
540 std::string message =
"Invalid line at line " + std::to_string(
lineCounter) +
541 ". Expected " + std::to_string(
fields.size()) +
542 " columns , but got " + std::to_string(results.size()) +
547 std::stringstream ss(line);
552 else if(it->name==
"Ek")
553 {
ReadValue(ss, Ek); Ek *= (CLHEP::GeV * it->unit);}
554 else if(it->name==
"E")
555 {
ReadValue(ss, E); E *= (CLHEP::GeV * it->unit);}
556 else if(it->name==
"P")
557 {
ReadValue(ss, P); P *= (CLHEP::GeV * it->unit);}
558 else if(it->name==
"t")
559 {
ReadValue(ss, t); t *= (CLHEP::s * it->unit); tdef =
true;}
560 else if(it->name==
"x")
561 {
ReadValue(ss, x); x *= (CLHEP::m * it->unit);}
562 else if(it->name==
"y")
563 {
ReadValue(ss, y); y *= (CLHEP::m * it->unit);}
564 else if(it->name==
"z")
565 {
ReadValue(ss, z); z *= (CLHEP::m * it->unit);}
566 else if(it->name==
"xp") {
ReadValue(ss, xp); xp *= ( CLHEP::radian * it->unit );}
567 else if(it->name==
"yp") {
ReadValue(ss, yp); yp *= ( CLHEP::radian * it->unit );}
568 else if(it->name==
"zp") {
ReadValue(ss, zp); zp *= ( CLHEP::radian * it->unit ); zpdef =
true;}
569 else if(it->name==
"pdgid")
572 updateParticleDefinition =
true;
574 else if (it->name ==
"S")
577 z *= CLHEP::m * it->unit;
579 else if(it->name==
"weight")
592 {updateParticleDefinition =
true;}
609 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
610 G4ParticleDefinition* particleDef =
nullptr;
613 {particleDef = particleTable->FindParticle(type);}
616 G4IonTable* ionTable = particleTable->GetIonTable();
617 G4int ionA, ionZ, ionLevel;
619 G4IonTable::GetNucleusByEncoding(type, ionZ, ionA, ionE, ionLevel);
625 {
throw BDSException(
"BDSBunchUserFile> Particle \"" + std::to_string(type) +
"\" not found");}
646template <
typename Type>
665 else if (fmt ==
"GeV")
667 else if (fmt ==
"MeV")
669 else if (fmt ==
"KeV")
671 else if (fmt ==
"eV")
674 {
throw BDSException(__METHOD_NAME__,
"Unrecognised energy unit! " +fmt);}
683 else if (fmt ==
"cm")
685 else if (fmt ==
"mm")
687 else if (fmt ==
"mum" || fmt ==
"um")
689 else if (fmt ==
"nm")
692 {
throw BDSException(__METHOD_NAME__,
"Unrecognised length unit! " + fmt);}
701 else if (fmt ==
"mrad")
703 else if (fmt ==
"murad" || fmt ==
"urad")
705 else if (fmt ==
"nrad")
708 {
throw BDSException(__METHOD_NAME__,
"Unrecognised angle unit! " + fmt);}
717 else if (fmt ==
"ms")
719 else if (fmt ==
"mus" || fmt ==
"us")
721 else if (fmt ==
"ns")
723 else if (fmt ==
"mm/c")
724 {
unit = (CLHEP::mm/CLHEP::c_light)/CLHEP::s;}
725 else if (fmt ==
"nm/c")
726 {
unit = (CLHEP::nm/CLHEP::c_light)/CLHEP::s;}
728 {
throw BDSException(__METHOD_NAME__,
"Unrecognised time unit! " + fmt);}
An intermediate layer for any bunch classes that are file based.
unsigned long long int nEventsInFile
The number of entries in the file loaded.
virtual void SetOptions(const BDSParticleDefinition *beamParticle, const GMAD::Beam &beam, const BDSBunchType &distrType, G4Transform3D beamlineTransformIn=G4Transform3D::Identity, const G4double beamlineS=0)
Pull out the relevant options and then pass through to BDSBunch::SetOptions().
A bunch distribution that reads a user specified column file.
void CheckConflictingParameters(const std::set< G4String > &s) const
Check conflicting columns aren't specified in file, e.g. P and Ek. Throw exception if wrong.
virtual void RecreateAdvanceToEvent(G4int eventOffset)
G4double ffact
Cache of flip factor from global constants.
void ReadValue(std::stringstream &stream, Type &value)
G4int lineCounter
Line counter.
G4long CountNLinesValidDataInFile()
void SkipNLinesSkip(G4bool usualPrintOut=true)
T InputBunchFile
The file handler. Templated as could be std::ifstream or igzstream for example.
G4bool SkippableLine(const std::string &line) const
Return true if a line is all whitespace or is commented out (starts with '#').
virtual BDSParticleCoordsFullGlobal GetNextParticleValid(G4int maxTries)
G4bool changingParticleType
Whether the particle type is a column.
void skip(std::stringstream &stream, G4int nvalues)
void SkipNLinesIgnoreIntoFile(G4bool usualPrintOut=true)
void ParseFileFormat()
Parse the column tokens and units factors.
G4long nlinesIgnore
Number of lines that will be ignored at the start the file.
virtual void SetOptions(const BDSParticleDefinition *beamParticle, const GMAD::Beam &beam, const BDSBunchType &distrType, G4Transform3D beamlineTransformIn=G4Transform3D::Identity, const G4double beamlineS=0)
Pull out the relevant options and then pass through to BDSBunch::SetOptions().
G4String distrFilePath
Bunch file including absolute path.
G4String bunchFormat
Format of the file.
G4long nlinesSkip
Number of lines that will be skipped after the nlinesIgnore.
G4double particleMass
Cache of nominal beam particle mass.
std::list< Doublet > fields
List of variables to parse on each line.
void CloseBunchFile()
Close the file handler.
virtual void CheckParameters()
virtual BDSParticleCoordsFull GetNextParticleLocal()
Get the next particle.
virtual void Initialise()
Open the file and skip lines.
G4String distrFile
Bunch file.
void OpenBunchFile()
Open the file and check it's open.
G4bool anEnergyCoordinateInUse
Whether Et, Ek or P are in the columns.
G4double Yp0
Centre of distributions.
virtual void RecreateAdvanceToEvent(G4int eventOffset)
G4bool particleDefinitionHasBeenUpdated
G4bool useCurvilinear
Whether to ignore z and use s and transform for curvilinear coordinates.
G4double Z0
Centre of distributions.
G4double X0
Centre of distributions.
G4double Xp0
Centre of distributions.
G4double E0
Centre of distributions.
G4double Y0
Centre of distributions.
BDSParticleCoordsFullGlobal GetNextParticle()
static G4double CalculateZp(G4double xp, G4double yp, G4double Zp0)
Calculate zp safely based on other components.
BDSParticleDefinition * particleDefinition
Particle definition for bunch - this class owns it.
virtual void CheckParameters()
General exception with possible name of object and message.
static BDSGlobalConstants * Instance()
Access method.
Class to parse an ion particle definition.
G4double ExcitationEnergy() const
Accessor.
A set of particle coordinates in both local and global.
A set of particle coordinates including energy and weight.
Wrapper for particle definition.
G4double Mass() const
Accessor.
void SetEnergies(G4double totalEnergyIn, G4double kineticEnergyIn, G4double momentumIn)
G4double TotalEnergy() const
Accessor.
Improve type-safety of native enum data type in C++.
std::string distrFileFormat
beam parameters
std::string distrFile
beam parameters
int nlinesIgnore
Ignore first lines in the input bunch file.
int nlinesSkip
Number of event lines to skip after the ignore lines.
bool distrFileMatchLength
beam parameters
Return either G4Tubs or G4CutTubs depending on flat face.
G4String GetFullPath(G4String filename, bool excludeNameFromPath=false, bool useCWDForPrefix=false)
G4double ParseTimeUnit(const G4String &fmt)
G4double ParseEnergyUnit(const G4String &fmt)
G4double ParseAngleUnit(const G4String &fmt)
G4bool IsFinite(G4double value, G4double tolerance=std::numeric_limits< double >::epsilon())
G4double ParseLengthUnit(const G4String &fmt)
Struct for name and unit pair.
G4double unit
relative to SI units, i.e. mm=0.001 etc.