GENIEGenerator
Loading...
Searching...
No Matches
HNLDecayUtils.cxx
Go to the documentation of this file.
1//____________________________________________________________________________
2/*
3 Copyright (c) 2003-2025, The GENIE Collaboration
4 For the full text of the license visit http://copyright.genie-mc.org
5
6 Author: John Plows <komninos-john.plows \at physics.ox.ac.uk>
7 University of Oxford
8
9 Costas Andreopoulos <c.andreopoulos \at cern.ch>
10 University of Liverpool
11*/
12//____________________________________________________________________________
13
14#include "Framework/Messenger/Messenger.h"
15#include "Framework/ParticleData/PDGCodes.h"
16#include "Framework/ParticleData/PDGUtils.h"
17#include "Framework/ParticleData/PDGLibrary.h"
18#include "Framework/Utils/PrintUtils.h"
19#include "Physics/BeamHNL/HNLDecayUtils.h"
20
21using namespace genie;
22using namespace genie::hnl;
23
24//____________________________________________________________________________
25string genie::utils::hnl::ProdAsString(HNLProd_t hnlprod)
26{
27 switch(hnlprod) {
28
29 case (kHNLProdPion2Muon):
30 return "pi+- --> N + mu+-";
31 break;
32
33 case (kHNLProdPion2Electron):
34 return "pi+- --> N + e+-";
35 break;
36
37 case (kHNLProdKaon2Muon):
38 return "K+- --> N + mu+-";
39 break;
40
41 case (kHNLProdKaon2Electron):
42 return "K+- --> N + e+-";
43 break;
44
45 case (kHNLProdKaon3Muon):
46 return "K+- --> N + mu+- + pi0";
47 break;
48
49 case (kHNLProdKaon3Electron):
50 return "K+- --> N + e+- + pi0";
51 break;
52
53 case (kHNLProdNeuk3Muon):
54 return "K0L --> N + mu+- + pi-+";
55 break;
56
57 case (kHNLProdNeuk3Electron):
58 return "K0L --> N + e+- + pi-+";
59 break;
60
61 case (kHNLProdMuon3Numu):
62 return "mu-+ --> N + numu(bar) + e-+";
63 break;
64
65 case (kHNLProdMuon3Nue):
66 return "mu-+ --> N + nue(bar) + e-+";
67 break;
68
69 case (kHNLProdMuon3Nutau):
70 return "mu-+ --> N + nutau(bar) + e-+";
71 break;
72
73 default:
74 return "Invalid HNL production mode!";
75 }
76}
77//____________________________________________________________________________
78string genie::utils::hnl::AsString(HNLDecayMode_t hnldm)
79{
80 switch(hnldm) {
81
82 case (kHNLDcyNull):
83 return "Unspecified";
84 break;
85
86 case (kHNLDcyPiMu):
87 return "N -> pi+- mu-+";
88 break;
89
90 case (kHNLDcyPiE):
91 return "N -> pi+- e-+";
92 break;
93
94 case (kHNLDcyPi0Nu):
95 return "N -> pi0 v";
96 break;
97
98 case (kHNLDcyNuNuNu):
99 return "N -> v v v";
100 break;
101
102 case (kHNLDcyNuMuMu):
103 return "N -> v mu+ mu-";
104 break;
105
106 case (kHNLDcyNuEE):
107 return "N -> v e+ e-";
108 break;
109
110 case (kHNLDcyNuMuE):
111 return "N -> v mu+- e-+";
112 break;
113
114 case (kHNLDcyPiPi0E):
115 return "N -> pi+- pi0 e-+";
116 break;
117
118 case (kHNLDcyPiPi0Mu):
119 return "N -> pi+- pi0 mu-+";
120 break;
121
122 case (kHNLDcyPi0Pi0Nu):
123 return "N -> v pi0 pi0";
124 break;
125
126 case (kHNLDcyTEST):
127 return "N -> v v";
128 break;
129
130 default:
131 return "Invalid HNL decay mode!";
132 }
133 //return "Invalid HNL decay mode!";
134}
135//____________________________________________________________________________
136bool genie::utils::hnl::IsProdKinematicallyAllowed(HNLProd_t hnlprod)
137{
138// Check the input mass of the HNL (M) against the sum of the masses of the
139// particles to be produced in the input production mode
140
141 PDGCodeList decay_products = ProductionProductList(hnlprod);
142
143 PDGLibrary * pdglib = PDGLibrary::Instance();
144
145 double Msum = 0.; double Mpar = 0.;
146 PDGCodeList::const_iterator it = decay_products.begin();
147 for ( ; it != decay_products.end(); ++it)
148 {
149 int pdg_code = *it;
150 TParticlePDG * p = pdglib->Find(pdg_code);
151 if(p) {
152 if( it == decay_products.begin() ) Mpar = p->Mass();
153 else Msum += p->Mass();
154 } else {
155 LOG("HNL", pERROR)
156 << "Decay list includes particle with unrecognised PDG code: "
157 << pdg_code;
158 }
159 }
160
161 return (Mpar > Msum);
162}
163//____________________________________________________________________________
164bool genie::utils::hnl::IsKinematicallyAllowed(HNLDecayMode_t hnldm, double M)
165{
166// Check the input mass of the HNL (M) against the sum of the masses of the
167// particles to be produced in the input decay
168
169 PDGCodeList decay_products = DecayProductList(hnldm);
170
171 PDGLibrary * pdglib = PDGLibrary::Instance();
172
173 double Msum = 0.;
174 PDGCodeList::const_iterator it = decay_products.begin();
175 for ( ; it != decay_products.end(); ++it)
176 {
177 int pdg_code = *it;
178 TParticlePDG * p = pdglib->Find(pdg_code);
179 if(p) {
180 Msum += p->Mass();
181 } else {
182 LOG("HNL", pERROR)
183 << "Decay list includes particle with unrecognised PDG code: "
184 << pdg_code;
185 }
186 }
187
188 return (M > Msum);
189}
190//____________________________________________________________________________
191PDGCodeList genie::utils::hnl::ProductionProductList(HNLProd_t hnldm)
192{
193 // 0th element is parent
194 // 1st is HNL
195 // the rest are the HNL's siblings
196 bool allow_duplicate = true;
197 PDGCodeList decay_products(allow_duplicate);
198
199 switch(hnldm) {
200
201 case(kHNLProdPion2Muon):
202 decay_products.push_back(kPdgPiP);
203 decay_products.push_back(kPdgHNL);
204 decay_products.push_back(kPdgMuon);
205 break;
206
207 case(kHNLProdPion2Electron):
208 decay_products.push_back(kPdgPiP);
209 decay_products.push_back(kPdgHNL);
210 decay_products.push_back(kPdgElectron);
211 break;
212
213 case(kHNLProdKaon2Muon):
214 decay_products.push_back(kPdgKP);
215 decay_products.push_back(kPdgHNL);
216 decay_products.push_back(kPdgMuon);
217 break;
218
219 case(kHNLProdKaon2Electron):
220 decay_products.push_back(kPdgKP);
221 decay_products.push_back(kPdgHNL);
222 decay_products.push_back(kPdgElectron);
223 break;
224
225 case(kHNLProdKaon3Muon):
226 decay_products.push_back(kPdgKP);
227 decay_products.push_back(kPdgHNL);
228 decay_products.push_back(kPdgMuon);
229 decay_products.push_back(kPdgPi0);
230 break;
231
232 case(kHNLProdKaon3Electron):
233 decay_products.push_back(kPdgKP);
234 decay_products.push_back(kPdgHNL);
235 decay_products.push_back(kPdgElectron);
236 decay_products.push_back(kPdgPi0);
237 break;
238
239 case(kHNLProdNeuk3Muon):
240 decay_products.push_back(kPdgK0L);
241 decay_products.push_back(kPdgHNL);
242 decay_products.push_back(kPdgMuon);
243 decay_products.push_back(kPdgPiP);
244 break;
245
246 case(kHNLProdNeuk3Electron):
247 decay_products.push_back(kPdgK0L);
248 decay_products.push_back(kPdgHNL);
249 decay_products.push_back(kPdgElectron);
250 decay_products.push_back(kPdgPiP);
251 break;
252
253 case(kHNLProdMuon3Numu):
254 decay_products.push_back(kPdgMuon);
255 decay_products.push_back(kPdgHNL);
256 decay_products.push_back(kPdgElectron);
257 decay_products.push_back(kPdgAntiNuMu);
258 break;
259
260 case(kHNLProdMuon3Nue):
261 decay_products.push_back(kPdgMuon);
262 decay_products.push_back(kPdgHNL);
263 decay_products.push_back(kPdgElectron);
264 decay_products.push_back(kPdgAntiNuE);
265 break;
266
267 case(kHNLProdMuon3Nutau):
268 decay_products.push_back(kPdgMuon);
269 decay_products.push_back(kPdgHNL);
270 decay_products.push_back(kPdgElectron);
271 decay_products.push_back(kPdgAntiNuTau);
272 break;
273
274 default :
275 break;
276 }
277
278 return decay_products;
279}
280//____________________________________________________________________________
281PDGCodeList genie::utils::hnl::DecayProductList(HNLDecayMode_t hnldm)
282{
283 bool allow_duplicate = true;
284 PDGCodeList decay_products(allow_duplicate);
285
286 switch(hnldm) {
287
288 case(kHNLDcyPiMu):
289 decay_products.push_back(kPdgPiP);
290 decay_products.push_back(kPdgMuon);
291 break;
292
293 case(kHNLDcyPiE):
294 decay_products.push_back(kPdgPiP);
295 decay_products.push_back(kPdgElectron);
296 break;
297
298 case(kHNLDcyPi0Nu):
299 decay_products.push_back(kPdgPi0);
300 decay_products.push_back(kPdgNuMu); // could be nue or nutau too!
301 break;
302
303 case(kHNLDcyNuNuNu):
304 decay_products.push_back(kPdgNuE);
305 decay_products.push_back(kPdgNuMu);
306 decay_products.push_back(kPdgNuTau); // again, any permutation of {e,mu,tau}^3 works
307 break;
308
309 case(kHNLDcyNuMuMu):
310 decay_products.push_back(kPdgNuMu);
311 decay_products.push_back(kPdgMuon);
312 decay_products.push_back(kPdgAntiMuon);
313 break;
314
315 case(kHNLDcyNuEE):
316 decay_products.push_back(kPdgNuE);
317 decay_products.push_back(kPdgElectron);
318 decay_products.push_back(kPdgPositron);
319 break;
320
321 case(kHNLDcyNuMuE):
322 decay_products.push_back(kPdgNuMu);
323 decay_products.push_back(kPdgMuon);
324 decay_products.push_back(kPdgPositron);
325 break;
326
327 case(kHNLDcyPiPi0E):
328 decay_products.push_back(kPdgPiP);
329 decay_products.push_back(kPdgPi0);
330 decay_products.push_back(kPdgElectron);
331 break;
332
333 case(kHNLDcyPiPi0Mu):
334 decay_products.push_back(kPdgPiP);
335 decay_products.push_back(kPdgPi0);
336 decay_products.push_back(kPdgMuon);
337 break;
338
339 case(kHNLDcyPi0Pi0Nu):
340 decay_products.push_back(kPdgPi0);
341 decay_products.push_back(kPdgPi0);
342 decay_products.push_back(kPdgNuMu);
343 break;
344
345 case(kHNLDcyTEST):
346 decay_products.push_back(kPdgNuMu);
347 decay_products.push_back(kPdgAntiNuMu);
348 break;
349
350 default :
351 break;
352 }
353
354 return decay_products;
355}
356//____________________________________________________________________________
357// see e.g. Physics/HadronTensors/TabulatedLabFrameHadronTensor.cxx for this code
358int genie::utils::hnl::GetCfgInt( string file_id, string set_name, string par_name )
359{
360 const genie::Registry * tmpReg = genie::AlgConfigPool::Instance()
361 ->CommonList( file_id, set_name );
362 int param = tmpReg->GetInt( par_name );
363 return param;
364}
365//____________________________________________________________________________
366// see e.g. Tools/Flux/GNuMIFlux.cxx for this code
367std::vector<int> genie::utils::hnl::GetCfgIntVec( string file_id, string set_name, string par_name )
368{
369 // get vector as string first
370 const genie::Registry * tmpReg = genie::AlgConfigPool::Instance()
371 ->CommonList( file_id, set_name );
372 string str = tmpReg->GetString( par_name );
373
374 // turn string into vector<int>
375 // be liberal about separators, users might punctuate for clarity
376 std::vector<string> strtokens = genie::utils::str::Split( str, " ,;:()[]=" );
377 std::vector<int> vect;
378 size_t ntok = strtokens.size();
379
380 for( size_t i=0; i < ntok; ++i ){
381 std::string trimmed = utils::str::TrimSpaces( strtokens[i] );
382 if( " " == trimmed || "" == trimmed ) continue; // skip empty strings
383 int val = strtod( trimmed.c_str(), (char **)NULL );
384 vect.push_back( val );
385 }
386
387 return vect;
388}
389//____________________________________________________________________________
390double genie::utils::hnl::GetCfgDouble( string file_id, string set_name, string par_name )
391{
392 const genie::Registry * tmpReg = genie::AlgConfigPool::Instance()
393 ->CommonList( file_id, set_name );
394 double param = tmpReg->GetDouble( par_name );
395 return param;
396}
397//____________________________________________________________________________
398std::vector<double> genie::utils::hnl::GetCfgDoubleVec( string file_id, string set_name, string par_name )
399{
400 // get vector as string first
401 const genie::Registry * tmpReg = genie::AlgConfigPool::Instance()
402 ->CommonList( file_id, set_name );
403 string str = tmpReg->GetString( par_name );
404
405 // turn string into vector<double>
406 // be liberal about separators, users might punctuate for clarity
407 std::vector<string> strtokens = genie::utils::str::Split( str, " ,;:()[]=" );
408 std::vector<double> vect;
409 size_t ntok = strtokens.size();
410
411 for( size_t i=0; i < ntok; ++i ){
412 std::string trimmed = utils::str::TrimSpaces( strtokens[i] );
413 if( " " == trimmed || "" == trimmed ) continue; // skip empty strings
414 double val = strtod( trimmed.c_str(), (char **)NULL );
415 vect.push_back( val );
416 }
417
418 return vect;
419}
420//____________________________________________________________________________
421bool genie::utils::hnl::GetCfgBool( string file_id, string set_name, string par_name )
422{
423 const genie::Registry * tmpReg = genie::AlgConfigPool::Instance()
424 ->CommonList( file_id, set_name );
425 bool param = tmpReg->GetBool( par_name );
426 return param;
427}
428//____________________________________________________________________________
429std::vector<bool> genie::utils::hnl::GetCfgBoolVec( string file_id, string set_name, string par_name )
430{
431 // get vector as string first
432 const genie::Registry * tmpReg = genie::AlgConfigPool::Instance()
433 ->CommonList( file_id, set_name );
434 string str = tmpReg->GetString( par_name );
435
436 // turn string into vector<bool>
437 // be liberal about separators, users might punctuate for clarity
438 std::vector<string> strtokens = genie::utils::str::Split( str, " ,;:()[]=" );
439 std::vector<bool> vect;
440 size_t ntok = strtokens.size();
441
442 for( size_t i=0; i < ntok; ++i ){
443 std::string trimmed = utils::str::TrimSpaces( strtokens[i] );
444 if( " " == trimmed || "" == trimmed ) continue; // skip empty strings
445 bool val = strtod( trimmed.c_str(), (char **)NULL );
446 vect.push_back( val );
447 }
448
449 return vect;
450}
451//____________________________________________________________________________
452string genie::utils::hnl::GetCfgString( string file_id, string set_name, string par_name )
453{
454 const genie::Registry * tmpReg = genie::AlgConfigPool::Instance()
455 ->CommonList( file_id, set_name );
456 string param = tmpReg->GetString( par_name );
457 return param;
458}
459//____________________________________________________________________________
pERROR
#define pERROR
Definition Messenger.h:59
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils.
genie::AlgConfigPool::Instance
static AlgConfigPool * Instance()
Definition AlgConfigPool.cxx:85
genie::AlgConfigPool::CommonList
Registry * CommonList(const string &file_id, const string &set_name) const
Definition AlgConfigPool.cxx:744
genie::PDGCodeList::push_back
void push_back(int pdg_code)
Definition PDGCodeList.cxx:58
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition PDGLibrary.cxx:68
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition PDGLibrary.cxx:86
genie::Registry::GetDouble
RgDbl GetDouble(RgKey key) const
Definition Registry.cxx:474
genie::Registry::GetString
RgStr GetString(RgKey key) const
Definition Registry.cxx:481
genie::Registry::GetInt
RgInt GetInt(RgKey key) const
Definition Registry.cxx:467
genie::Registry::GetBool
RgBool GetBool(RgKey key) const
Definition Registry.cxx:460
genie::hnl::HNLProd_t
enum genie::hnl::t_HNLProd HNLProd_t
genie::hnl::HNLDecayMode_t
enum genie::hnl::EHNLDecayMode HNLDecayMode_t
genie::utils::hnl::DecayProductList
PDGCodeList DecayProductList(genie::hnl::HNLDecayMode_t hnldm)
genie::utils::hnl::GetCfgBoolVec
std::vector< bool > GetCfgBoolVec(string file_id, string set_name, string par_name)
genie::utils::hnl::GetCfgDoubleVec
std::vector< double > GetCfgDoubleVec(string file_id, string set_name, string par_name)
genie::utils::hnl::GetCfgString
std::string GetCfgString(string file_id, string set_name, string par_name)
genie::utils::hnl::GetCfgInt
int GetCfgInt(string file_id, string set_name, string par_name)
genie::utils::hnl::ProductionProductList
PDGCodeList ProductionProductList(genie::hnl::HNLProd_t hnldm)
genie::utils::hnl::GetCfgBool
bool GetCfgBool(string file_id, string set_name, string par_name)
genie::utils::hnl::GetCfgDouble
double GetCfgDouble(string file_id, string set_name, string par_name)
genie::utils::hnl::ProdAsString
string ProdAsString(genie::hnl::HNLProd_t hnlprod)
genie::utils::hnl::AsString
string AsString(genie::hnl::HNLDecayMode_t hnldm)
genie::utils::hnl::IsKinematicallyAllowed
bool IsKinematicallyAllowed(genie::hnl::HNLDecayMode_t hnldm, double Mhnl)
genie::utils::hnl::IsProdKinematicallyAllowed
bool IsProdKinematicallyAllowed(genie::hnl::HNLProd_t hnlprod)
genie::utils::hnl::GetCfgIntVec
std::vector< int > GetCfgIntVec(string file_id, string set_name, string par_name)
genie::utils::str
Utilities for string manipulation.
genie::utils::str::TrimSpaces
string TrimSpaces(string input)
Definition StringUtils.cxx:18
genie::utils::str::Split
vector< string > Split(string input, string delim)
Definition StringUtils.cxx:36
genie
THE MAIN GENIE PROJECT NAMESPACE
Definition AlgCmp.h:25
genie::kPdgAntiMuon
const int kPdgAntiMuon
Definition PDGCodes.h:38
genie::kPdgPi0
const int kPdgPi0
Definition PDGCodes.h:160
genie::kPdgK0L
const int kPdgK0L
Definition PDGCodes.h:176
genie::kPdgKP
const int kPdgKP
Definition PDGCodes.h:172
genie::kPdgAntiNuE
const int kPdgAntiNuE
Definition PDGCodes.h:29
genie::kPdgAntiNuTau
const int kPdgAntiNuTau
Definition PDGCodes.h:33
genie::kPdgNuE
const int kPdgNuE
Definition PDGCodes.h:28
genie::kPdgNuTau
const int kPdgNuTau
Definition PDGCodes.h:32
genie::kPdgHNL
const int kPdgHNL
Definition PDGCodes.h:224
genie::kPdgAntiNuMu
const int kPdgAntiNuMu
Definition PDGCodes.h:31
genie::kPdgPiP
const int kPdgPiP
Definition PDGCodes.h:158
genie::kPdgMuon
const int kPdgMuon
Definition PDGCodes.h:37
genie::kPdgNuMu
const int kPdgNuMu
Definition PDGCodes.h:30
genie::kPdgPositron
const int kPdgPositron
Definition PDGCodes.h:36
genie::kPdgElectron
const int kPdgElectron
Definition PDGCodes.h:35
Generated on for GENIEGenerator by doxygen 1.14.0