GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HNLDecayUtils.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, 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 
21 using namespace genie;
22 using namespace genie::hnl;
23 
24 //____________________________________________________________________________
25 string 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 //____________________________________________________________________________
78 string 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 //____________________________________________________________________________
136 bool 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 //____________________________________________________________________________
164 bool 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 //____________________________________________________________________________
191 PDGCodeList 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 //____________________________________________________________________________
281 PDGCodeList 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
358 int 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
367 std::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 //____________________________________________________________________________
390 double 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 //____________________________________________________________________________
398 std::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 //____________________________________________________________________________
421 bool 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 //____________________________________________________________________________
429 std::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 //____________________________________________________________________________
452 string 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 //____________________________________________________________________________
genie::hnl::HNLProd_t
enum genie::hnl::t_HNLProd HNLProd_t
genie::utils::hnl::GetCfgDouble
double GetCfgDouble(string file_id, string set_name, string par_name)
genie::kPdgNuE
const int kPdgNuE
Definition: PDGCodes.h:28
pERROR
#define pERROR
Definition: Messenger.h:59
genie::utils::hnl::AsString
string AsString(genie::hnl::HNLDecayMode_t hnldm)
genie::utils::hnl::ProductionProductList
PDGCodeList ProductionProductList(genie::hnl::HNLProd_t hnldm)
genie::kPdgAntiNuE
const int kPdgAntiNuE
Definition: PDGCodes.h:29
genie::kPdgNuMu
const int kPdgNuMu
Definition: PDGCodes.h:30
genie::kPdgAntiMuon
const int kPdgAntiMuon
Definition: PDGCodes.h:38
genie::Registry::GetDouble
RgDbl GetDouble(RgKey key) const
Definition: Registry.cxx:474
genie::kPdgElectron
const int kPdgElectron
Definition: PDGCodes.h:35
genie::utils::hnl::GetCfgString
std::string GetCfgString(string file_id, string set_name, string par_name)
genie::Registry::GetInt
RgInt GetInt(RgKey key) const
Definition: Registry.cxx:467
genie::utils::hnl::GetCfgBoolVec
std::vector< bool > GetCfgBoolVec(string file_id, string set_name, string par_name)
genie::AlgConfigPool::CommonList
Registry * CommonList(const string &file_id, const string &set_name) const
Definition: AlgConfigPool.cxx:744
genie::PDGCodeList
A list of PDG codes.
Definition: PDGCodeList.h:32
genie::utils::hnl::GetCfgIntVec
std::vector< int > GetCfgIntVec(string file_id, string set_name, string par_name)
genie::utils::hnl::DecayProductList
PDGCodeList DecayProductList(genie::hnl::HNLDecayMode_t hnldm)
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::kPdgHNL
const int kPdgHNL
Definition: PDGCodes.h:224
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::hnl::HNLDecayMode_t
enum genie::hnl::EHNLDecayMode HNLDecayMode_t
genie::utils::hnl::ProdAsString
string ProdAsString(genie::hnl::HNLProd_t hnlprod)
genie::kPdgK0L
const int kPdgK0L
Definition: PDGCodes.h:176
genie::kPdgAntiNuTau
const int kPdgAntiNuTau
Definition: PDGCodes.h:33
genie::kPdgAntiNuMu
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
genie::utils::str::TrimSpaces
string TrimSpaces(string input)
Definition: StringUtils.cxx:18
genie::utils::hnl::GetCfgDoubleVec
std::vector< double > GetCfgDoubleVec(string file_id, string set_name, string par_name)
genie::utils::hnl::GetCfgBool
bool GetCfgBool(string file_id, string set_name, string par_name)
genie::kPdgNuTau
const int kPdgNuTau
Definition: PDGCodes.h:32
genie::Registry::GetString
RgStr GetString(RgKey key) const
Definition: Registry.cxx:481
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::utils::hnl::IsProdKinematicallyAllowed
bool IsProdKinematicallyAllowed(genie::hnl::HNLProd_t hnlprod)
genie::utils::str::Split
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
genie::PDGLibrary
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
genie::Registry
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
genie::Registry::GetBool
RgBool GetBool(RgKey key) const
Definition: Registry.cxx:460
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
genie::kPdgMuon
const int kPdgMuon
Definition: PDGCodes.h:37
genie::kPdgPositron
const int kPdgPositron
Definition: PDGCodes.h:36
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
genie::utils::hnl::IsKinematicallyAllowed
bool IsKinematicallyAllowed(genie::hnl::HNLDecayMode_t hnldm, double Mhnl)
genie::PDGCodeList::push_back
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
genie::AlgConfigPool::Instance
static AlgConfigPool * Instance()
Definition: AlgConfigPool.cxx:85
genie::utils::hnl::GetCfgInt
int GetCfgInt(string file_id, string set_name, string par_name)
Generated on Wed Jul 17 2024 14:46:42 for GENIEGenerator by   doxygen 1.8.5