GENIEGenerator
Loading...
Searching...
No Matches
Public Member Functions | Private Member Functions | Private Attributes | List of all members
genie::hnl::Decayer Class Reference

Heavy Neutral Lepton final-state product generator. More...

#include <HNLDecayer.h>

Inheritance diagram for genie::hnl::Decayer:
[legend]
Collaboration diagram for genie::hnl::Decayer:
[legend]

Public Member Functions

 Decayer ()
 Decayer (string config)
 ~Decayer ()
void ProcessEventRecord (GHepRecord *event) const
void Configure (const Registry &config)
void Configure (string config)
double GetHNLLifetime () const
double GetHNLMass () const
std::vector< double > GetHNLCouplings () const
bool IsHNLMajorana () const
std::string GetHNLInterestingChannels () const
std::vector< double > GetPGunOrigin () const
std::vector< double > GetPGunDOrigin () const
double GetPGunEnergy () const
std::vector< double > GetPGunDirection () const
std::vector< double > GetPGunDeviation () const
Public Member Functions inherited from genie::hnl::DecayRecordVisitorI
virtual ~DecayRecordVisitorI ()
Public Member Functions inherited from genie::EventRecordVisitorI
virtual ~EventRecordVisitorI ()
Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
virtual void FindConfig (void)
virtual const Registry & GetConfig (void) const
Registry * GetOwnedConfig (void)
virtual const AlgId & Id (void) const
 Get algorithm ID.
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status.
virtual bool AllowReconfig (void) const
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm.
virtual void SetId (const AlgId &id)
 Set algorithm ID.
virtual void SetId (string name, string config)
const Algorithm * SubAlg (const RgKey &registry_key) const
void AdoptConfig (void)
void AdoptSubstructure (void)
virtual void Print (ostream &stream) const
 Print algorithm info.

Private Member Functions

void LoadConfig (void)
void AddInitialState (GHepRecord *event) const
void GenerateDecayProducts (GHepRecord *event) const
void UpdateEventRecord (GHepRecord *event) const
void SetHNLCouplings (double Ue42, double Um42, double Ut42) const
void SetBeam2User (std::vector< double > translation, std::vector< double > rotation) const
void SetProdVtxPosition (const TLorentzVector &v4) const
void ReadCreationInfo (GHepRecord *event) const
genie::hnl::SimpleHNL GetHNLInstance () const
std::vector< double > * GenerateDecayPosition (GHepRecord *event) const
std::vector< double > * GenerateMomentum (GHepRecord *event) const
double CalcPolMag (int parPdg, int lepPdg, double M) const
double CalcPolMod (double polMag, int lepPdg, int hadPdg, double M) const
bool UnpolarisedDecay (TGenPhaseSpace &fPSG, PDGCodeList pdgv, double wm) const
bool PolarisedDecay (TGenPhaseSpace &fPSG, PDGCodeList pdgv, double wm, TVector3 vPolDir) const

Private Attributes

int fCurrInitStatePdg
genie::hnl::HNLDecayMode_t fCurrDecayMode
int fParentPdg
int fProdLepPdg
int fDecLepPdg
int fDecHadPdg
std::vector< double > fPolDir
bool fIsConfigLoaded = false
double fMass
double fUe42 = -1.0
double fUm42 = -1.0
double fUt42 = -1.0
bool fIsMajorana = false
bool fDoPol = false
std::vector< genie::hnl::HNLDecayMode_t > fIntChannels
int fChanBits [10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
std::vector< double > fB2UTranslation
double fTx = -1.0
double fTy = -1.0
double fTz = -1.0
std::vector< double > fB2URotation
double fR1 = -1.0
double fR2 = -1.0
double fR3 = -1.0
TH3D * fProdVtxHist = 0
TLorentzVector * fProdVtx = 0
TLorentzVector * fISMom = 0
double fPGOx = 0.0
double fPGOy = 0.0
double fPGOz = 0.0
double fPGDx = 0.0
double fPGDy = 0.0
double fPGDz = 0.0
double fPGE = 0.0
double fPGCx = 0.0
double fPGCy = 0.0
double fPGCz = 0.0
double fPGDTheta = 0.0
double fPGDPhi = 0.0
bool fGetCMFrameInstead = false

Additional Inherited Members

Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
static string BuildParamVectSizeKey (const std::string &comm_name)
static string BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
static string BuildParamMatRowSizeKey (const std::string &comm_name)
static string BuildParamMatColSizeKey (const std::string &comm_name)
Protected Member Functions inherited from genie::hnl::DecayRecordVisitorI
 DecayRecordVisitorI ()
 DecayRecordVisitorI (string name)
 DecayRecordVisitorI (string name, string config)
Protected Member Functions inherited from genie::EventRecordVisitorI
 EventRecordVisitorI ()
 EventRecordVisitorI (string name)
 EventRecordVisitorI (string name, string config)
Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 Algorithm (string name)
 Algorithm (string name, string config)
void Initialize (void)
void DeleteConfig (void)
void DeleteSubstructure (void)
Registry * ExtractLocalConfig (const Registry &in) const
Registry * ExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key.
template<class T>
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
template<class T>
bool GetParamDef (const RgKey &name, T &p, const T &def) const
template<class T>
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters.
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
template<class T>
int GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
 Handle to load matrix of parameters.
template<class T>
int GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
int GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership
int MergeTopRegistry (const Registry &r)
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships.
Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...)
AlgId fID
 algorithm name and configuration set
vector< Registry * > fConfVect
vector< bool > fOwnerships
 ownership for every registry in fConfVect
AlgStatus_t fStatus
 algorithm execution status
AlgMap * fOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool)

Detailed Description

Heavy Neutral Lepton final-state product generator.

Author
Costas Andreopoulos <c.andreopoulos \at cern.ch> University of Liverpool John Plows <komninos-john.plows \at physics.ox.ac.uk>
Created:\n February 10, 2020
License:\n Copyright (c) 2003-2025, The GENIE Collaboration
For the full text of the license visit http://copyright.genie-mc.org

Definition at line 41 of file HNLDecayer.h.

Constructor & Destructor Documentation

◆ Decayer() [1/2]

Decayer::Decayer ( )

Definition at line 37 of file HNLDecayer.cxx.

37 :
38DecayRecordVisitorI("genie::hnl::Decayer")
39{
40
41}

References genie::hnl::DecayRecordVisitorI::DecayRecordVisitorI().

◆ Decayer() [2/2]

Decayer::Decayer ( string config)

Definition at line 43 of file HNLDecayer.cxx.

43 :
44DecayRecordVisitorI("genie::hnl::Decayer",config)
45{
46
47}

References genie::hnl::DecayRecordVisitorI::DecayRecordVisitorI().

◆ ~Decayer()

Decayer::~Decayer ( )

Definition at line 49 of file HNLDecayer.cxx.

50{
51
52}

Member Function Documentation

◆ AddInitialState()

void Decayer::AddInitialState ( GHepRecord * event) const
private

Definition at line 73 of file HNLDecayer.cxx.

74{
75 std::vector< double > * prodVtx = 0;
76
77 Interaction * interaction = event->Summary();
78 InitialState * init_state = interaction->InitStatePtr();
79
80 TLorentzVector p4;
81 if( event->Particle(0) ){
82 // p4 was already set using HNLFluxCreator. No action needed.
83 // Read event vertex == HNL production vertex. We will find the decay vertex later.
84 p4 = *( init_state->GetProbeP4() );
85
86 prodVtx = new std::vector< double >();
87 prodVtx->emplace_back( event->Vertex()->X() );
88 prodVtx->emplace_back( event->Vertex()->Y() );
89 prodVtx->emplace_back( event->Vertex()->Z() );
90 prodVtx->emplace_back( event->Vertex()->T() );
91 } else {
92 std::vector< double > * p3HNL = this->GenerateMomentum( event );
93
94 double px = p3HNL->at(0);
95 double py = p3HNL->at(1);
96 double pz = p3HNL->at(2);
97 double E = interaction->InitState().ProbeE(kRfLab);
98
99 p4 = TLorentzVector( px, py, pz, E );
100
101 if( !event->Vertex() || (event->Vertex()->Vect()).Mag() == 0.0 )
102 prodVtx = this->GenerateDecayPosition( event );
103 else{
104 prodVtx = new std::vector< double >();
105 prodVtx->emplace_back( event->Vertex()->X() );
106 prodVtx->emplace_back( event->Vertex()->Y() );
107 prodVtx->emplace_back( event->Vertex()->Z() );
108 prodVtx->emplace_back( event->Vertex()->T() );
109 }
110 }
111
112 TLorentzVector v4( prodVtx->at(0), prodVtx->at(1), prodVtx->at(2), prodVtx->at(3) );
113
114 init_state->SetProbeP4( p4 );
115
116 int hpdg = interaction->InitState().ProbePdg();
117 if( !event->Particle(0) )
118 event->AddParticle(hpdg, kIStInitialState, 0,-1,-1,-1, p4, v4);
119}
genie::GHepRecord::Vertex
virtual TLorentzVector * Vertex(void) const
Definition GHepRecord.h:140
genie::GHepRecord::Particle
virtual GHepParticle * Particle(int position) const
Definition GHepRecord.cxx:104
genie::InitialState::GetProbeP4
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
Definition InitialState.cxx:325
genie::InitialState::SetProbeP4
void SetProbeP4(const TLorentzVector &P4)
Definition InitialState.cxx:183
genie::InitialState::ProbePdg
int ProbePdg(void) const
Definition InitialState.h:64
genie::InitialState::ProbeE
double ProbeE(RefFrame_t rf) const
Definition InitialState.cxx:384
genie::Interaction::InitStatePtr
InitialState * InitStatePtr(void) const
Definition Interaction.h:74
genie::Interaction::InitState
const InitialState & InitState(void) const
Definition Interaction.h:69
genie::hnl::Decayer::GenerateMomentum
std::vector< double > * GenerateMomentum(GHepRecord *event) const
Definition HNLDecayer.cxx:351
genie::hnl::Decayer::GenerateDecayPosition
std::vector< double > * GenerateDecayPosition(GHepRecord *event) const
Definition HNLDecayer.cxx:326
genie::kIStInitialState
@ kIStInitialState
Definition GHepStatus.h:29
genie::kRfLab
@ kRfLab
Definition RefFrame.h:26

References GenerateDecayPosition(), GenerateMomentum(), genie::InitialState::GetProbeP4(), genie::Interaction::InitState(), genie::Interaction::InitStatePtr(), genie::kIStInitialState, genie::kRfLab, genie::GHepRecord::Particle(), genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), genie::InitialState::SetProbeP4(), and genie::GHepRecord::Vertex().

Referenced by ProcessEventRecord().

◆ CalcPolMag()

double Decayer::CalcPolMag ( int parPdg,
int lepPdg,
double M ) const
private

Definition at line 671 of file HNLDecayer.cxx.

672{
673 PDGLibrary * pdgl = PDGLibrary::Instance();
674 double mPar = pdgl->Find( std::abs( parPdg ) )->Mass();
675 double mLep = pdgl->Find( std::abs( lepPdg ) )->Mass();
676
677 double num1 = mLep * mLep - M * M;
678 double num2 = TMath::Sqrt(utils::hnl::Kallen( mPar*mPar, M*M, mLep*mLep ));
679
680 double den1 = mPar*mPar*( mLep*mLep + M*M );
681 double den2 = mLep*mLep - M*M;
682
683 // pMag is a modulus, not a magnitude... not positive semi-definite. See Fig.4 in 1805.06419[hep-ph]
684 double pMag = -1.0 * num1*num2 / ( den1 - den2*den2 );
685 return pMag;
686}
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::utils::hnl::Kallen
double Kallen(double x, double y, double z)
Definition HNLKinUtils.h:37

References genie::PDGLibrary::Find(), genie::PDGLibrary::Instance(), and genie::utils::hnl::Kallen().

Referenced by PolarisedDecay().

◆ CalcPolMod()

double Decayer::CalcPolMod ( double polMag,
int lepPdg,
int hadPdg,
double M ) const
private

Definition at line 688 of file HNLDecayer.cxx.

689{
690 PDGLibrary * pdgl = PDGLibrary::Instance();
691 double mLep = pdgl->Find( std::abs( lepPdg ) )->Mass();
692 double mHad = pdgl->Find( std::abs( hadPdg ) )->Mass();
693
694 double num1 = M*M - mLep*mLep;
695 double num2 = TMath::Sqrt(utils::hnl::Kallen( M*M, mLep*mLep, mHad*mHad ));
696 double num3 = polMag;
697
698 double den1 = M*M - mLep*mLep;
699 double den2 = mHad*mHad*( M*M + mLep*mLep );
700
701 double pMod = num1*num2*num3 / ( den1*den1 - den2 );
702 return pMod;
703}

References genie::PDGLibrary::Find(), genie::PDGLibrary::Instance(), and genie::utils::hnl::Kallen().

Referenced by PolarisedDecay().

◆ Configure() [1/2]

void Decayer::Configure ( const Registry & config)
virtual

Configure the algorithm with an external registry The registry is merged with the top level registry if it is owned, Otherwise a copy of it is added with the highest priority

Reimplemented from genie::Algorithm.

Definition at line 407 of file HNLDecayer.cxx.

408{
409 Algorithm::Configure(config);
410 this->LoadConfig();
411}
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition Algorithm.cxx:62

References genie::Algorithm::Configure(), and LoadConfig().

◆ Configure() [2/2]

void Decayer::Configure ( string config)
virtual

Configure the algorithm from the AlgoConfigPool based on param_set string given in input An algorithm contains a vector of registries coming from different xml configuration files, which are loaded according a very precise prioriy This methods will load a number registries in order of priority: 1) "Tunable" parameter set from CommonParametes. This is loaded with the highest prioriry and it is designed to be used for tuning procedure Usage not expected from the user. 2) For every string defined in "CommonParame" the corresponding parameter set will be loaded from CommonParameter.xml 3) parameter set specified by the config string and defined in the xml file of the algorithm 4) if config is not "Default" also the Default parameter set from the same xml file will be loaded Effectively this avoids the repetion of a parameter when it is not changed in the requested configuration

Reimplemented from genie::Algorithm.

Definition at line 413 of file HNLDecayer.cxx.

414{
415 Algorithm::Configure(config);
416 this->LoadConfig();
417}

References genie::Algorithm::Configure(), and LoadConfig().

◆ GenerateDecayPosition()

std::vector< double > * Decayer::GenerateDecayPosition ( GHepRecord * event) const
private

Definition at line 326 of file HNLDecayer.cxx.

327{
328 // let's query *where* the HNL decayed from.
329 LOG( "HNL", pWARN )
330 << "\nYou are seeing this message because the input dk2nu files did not give position information (for some reason)..."
331 << "\nDistributing the production vertex at some point in a 1m-side box. This is not good, and results will likely be unphysical.";
332 fProdVtxHist = new TH3D( "dummy", "dummy", 100, -0.5, 0.5, 100, -0.5, 0.5, 100, -0.5, 0.5 );
333 assert( fProdVtxHist );
334
335 RandomGen * rnd = RandomGen::Instance();
336 double pvx = (rnd->RndGen()).Uniform( -0.5, 0.5 );
337 double pvy = (rnd->RndGen()).Uniform( -0.5, 0.5 );
338 double pvz = (rnd->RndGen()).Uniform( -0.5, 0.5 );
339 std::vector< double > * prodVtx = new std::vector< double >();
340 prodVtx->emplace_back( pvx ); prodVtx->emplace_back( pvy ); prodVtx->emplace_back( pvz );
341 LOG( "HNL", pDEBUG )
342 << "Production vertex at: ( " << prodVtx->at(0) << ", " << prodVtx->at(1) << ", " << prodVtx->at(2) << ") [cm]";
343
344 TLorentzVector v4( prodVtx->at(0), prodVtx->at(1), prodVtx->at(2), 0.0 );
345
346 SetProdVtxPosition( v4 );
347
348 return prodVtx;
349}
pDEBUG
#define pDEBUG
Definition Messenger.h:63
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
pWARN
#define pWARN
Definition Messenger.h:60
genie::RandomGen::Instance
static RandomGen * Instance()
Access instance.
Definition RandomGen.cxx:74
genie::RandomGen::RndGen
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition RandomGen.h:80
genie::hnl::Decayer::SetProdVtxPosition
void SetProdVtxPosition(const TLorentzVector &v4) const
Definition HNLDecayer.cxx:511

References fProdVtxHist, genie::RandomGen::Instance(), LOG, pDEBUG, pWARN, genie::RandomGen::RndGen(), and SetProdVtxPosition().

Referenced by AddInitialState().

◆ GenerateDecayProducts()

void Decayer::GenerateDecayProducts ( GHepRecord * event) const
private

Definition at line 121 of file HNLDecayer.cxx.

122{
123 LOG("HNL", pINFO) << "Generating decay...";
124 fDecLepPdg = 0; fDecHadPdg = 0;
125
126 // do we have nubar?
127 PDGCodeList pdgv0 = utils::hnl::DecayProductList(fCurrDecayMode);
128 int typeMod = ( fCurrInitStatePdg >= 0 ) ? 1 : -1;
129 if( fCurrDecayMode == kHNLDcyPi0Nu ){
130 fDecHadPdg = typeMod * kPdgPi0;
131 } else if( fCurrDecayMode != kHNLDcyNuNuNu ){
132 fDecHadPdg = typeMod * kPdgPiP;
133 }
134
135 if( event->Particle(0) && !fIsMajorana ){
136 typeMod = ( event->Particle(0)->Pdg() >= 0 ) ? 1 : -1;
137 } else if( event->Particle(0) && fIsMajorana ){
138 // equal probability to decay to 1 of 2 charge-conjugated states
139 RandomGen * Rng = RandomGen::Instance();
140 double rnd = Rng->RndGen().Uniform(0.0, 1.0);
141 typeMod = ( rnd >= 0.5 ) ? 1.0 : -1.0;
142 }
143 PDGCodeList pdgv(true);
144 for( std::vector<int>::iterator it = pdgv0.begin(); it != pdgv0.end(); ++it ){
145 int pdgc = *it;
146 int newpdgc = ( pdgc == genie::kPdgPi0 ) ? pdgc : typeMod * pdgc; // pi-0 is its own antiparticle
147 pdgv.push_back( newpdgc );
148 }
149
150 assert ( pdgv.size() > 1);
151
152 // if user wants to include polarisation effects, start prep now
153 /*
154 double fPolDirMag = 0.0;
155 if( fPolDir.size() == 3 ){
156 fPolDirMag = std::sqrt( ( fPolDir.at(0) * fPolDir.at(0) ) +
157 ( fPolDir.at(1) * fPolDir.at(1) ) +
158 ( fPolDir.at(2) * fPolDir.at(2) ) );
159 }
160 */
161 bool doPol = fDoPol && (fPolDir.size() == 3);
162
163 std::ostringstream asts;
164 if( !doPol ) asts << "Performing a phase space decay...";
165 else asts << "Performing a polarised decay with polarisation vector:"
166 << " ( " << fPolDir.at(0) << ", " << fPolDir.at(1) << ", " << fPolDir.at(2) << " )";
167 LOG("HNL", pINFO) << asts.str();
168
169 // Get the decay product masses
170
171 vector<int>::const_iterator pdg_iter;
172 int i = 0;
173 double * mass = new double[pdgv.size()];
174 double sum = 0;
175 for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
176 int pdgc = *pdg_iter;
177 double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
178 mass[i++] = m;
179 sum += m;
180 }
181
182 LOG("HNL", pINFO)
183 << "Decaying N = " << pdgv.size() << " particles / total mass = " << sum;
184
185 int hnl_id = 0;
186 GHepParticle * hnl = event->Particle(hnl_id);
187 assert(hnl);
188 TLorentzVector * p4d = hnl->GetP4(); TVector3 HNLBVec = p4d->BoostVector();
189 TLorentzVector * p4d_rest = (TLorentzVector *) p4d->Clone(); p4d_rest->Boost( -HNLBVec );
190 TLorentzVector * v4d = hnl->GetX4();
191
192 LOG("HNL", pINFO)
193 << "\nDecaying system p4 = " << utils::print::P4AsString(p4d)
194 << "\nin rest frame, p4 = " << utils::print::P4AsString(p4d_rest);
195
196 // Set the decay
197 TGenPhaseSpace fPhaseSpaceGenerator;
198 //bool permitted = fPhaseSpaceGenerator.SetDecay(*p4d, pdgv.size(), mass);
199 bool permitted = fPhaseSpaceGenerator.SetDecay(*p4d_rest, pdgv.size(), mass);
200 if(!permitted) {
201 LOG("HNL", pERROR)
202 << " *** Phase space decay is not permitted \n"
203 << " Total particle mass = " << sum << "\n"
204 << " Decaying system p4 = " << utils::print::P4AsString(p4d);
205 // clean-up
206 delete [] mass;
207 delete p4d; delete p4d_rest;
208 delete v4d;
209 // throw exception
210 genie::exceptions::EVGThreadException exception;
211 exception.SetReason("Decay not permitted kinematically");
212 exception.SwitchOnFastForward();
213 throw exception;
214 }
215
216 // Get the maximum weight
217 //double wmax = fPhaseSpaceGenerator.GetWtMax();
218 double wmax = -1;
219 for(int idec=0; idec<200; idec++) {
220 double w = fPhaseSpaceGenerator.Generate();
221 wmax = TMath::Max(wmax,w);
222 }
223 assert(wmax>0);
224 wmax *= 2;
225
226 LOG("HNL", pNOTICE)
227 << "Max phase space gen. weight @ current hadronic system: " << wmax;
228
229 // Generate a decay
230 bool decayFailed = false;
231 if( doPol && fCurrDecayMode != kHNLDcyNuNuNu ){
232 // if polarisation is on we must calculate the polarisation modulus for comparison
233 // for now, assume 2-body production and 2-body decay describes the pol modulus
234 // see arXiv:1805.06419[hep-ph]
235 TVector3 vecPolDir( fPolDir.at(0), fPolDir.at(1), fPolDir.at(2) );
236 LOG( "HNL", pDEBUG ) << "Doing a polarised decay...";
237 decayFailed = this->PolarisedDecay( fPhaseSpaceGenerator, pdgv, wmax, vecPolDir );
238 } else {
239 if( doPol && fCurrDecayMode == kHNLDcyNuNuNu ){
240 // no charged lepton here... warn the user about it, though
241 LOG( "HNL", pWARN )
242 << "Polarisation for invisible FS not implemented yet, defaulting to phase-space decay...";
243 }
244
245 LOG( "HNL", pDEBUG ) << "Doing a phase-space decay...";
246 decayFailed = this->UnpolarisedDecay( fPhaseSpaceGenerator, pdgv, wmax );
247 }
248 if( decayFailed ){
249 // clean up
250 delete [] mass;
251 delete p4d;
252 delete v4d;
253 // throw exception
254 genie::exceptions::EVGThreadException exception;
255 exception.SetReason("Couldn't select decay after N attempts");
256 exception.SwitchOnFastForward();
257 throw exception;
258 }
259
260 // Insert final state products into a TClonesArray of GHepParticle's
261 TLorentzVector v4(*v4d);
262 int idp = 0; int npip = 0, npi0 = 0, npim = 0;
263 for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
264 int pdgc = *pdg_iter;
265 TLorentzVector * p4fin = fPhaseSpaceGenerator.GetDecay(idp);
266 if( !fGetCMFrameInstead )
267 p4fin->Boost( HNLBVec ); // from rest to lab again
268 GHepStatus_t ist = kIStStableFinalState;
269 event->AddParticle(pdgc, ist, hnl_id,-1,-1,-1, *p4fin, v4);
270 switch( pdgc ){
271 case kPdgPiP: npip++; break;
272 case kPdgPi0: npi0++; break;
273 case kPdgPiM: npim++; break;
274 }
275 idp++;
276 }
277 Interaction * interaction = event->Summary();
278 interaction->ExclTagPtr()->SetNPions( npip, npi0, npim );
279
280 // Manually set up some mother-daughter links
281 // find primary (=leading) *charged!* lepton
282 double Elead = -1.0; int ilead = -1;
283 std::vector< int > lpdgv = { 11, 13, 15 };
284 for( std::vector<int>::iterator lit = lpdgv.begin(); lit != lpdgv.end(); ++lit ) {
285 GHepParticle * tmpPart = event->FindParticle( (*lit), kIStStableFinalState, 1 );
286 if( !tmpPart ) tmpPart = event->FindParticle( -1 * (*lit), kIStStableFinalState, 1 ); //antiparticle?
287 if( tmpPart ){
288 double tmpE = tmpPart->E();
289 if( tmpE > Elead ){ Elead = tmpE; ilead = event->ParticlePosition( tmpPart ); }
290 }
291 }
292 event->Particle( 0 )->SetFirstDaughter( ilead );
293 event->Particle( 0 )->SetFirstMother(-1); // why do I need to do this explicitly?
294 event->Particle( 0 )->SetLastMother(-1);
295
296 assert( event->Probe() );
297 if( !event->FinalStatePrimaryLepton() ){ // no charged lepton means invisible or pi0 nu or test
298 LOG( "HNL", pWARN )
299 << "No final state primary lepton for this event.";
300 assert( fCurrDecayMode == kHNLDcyPi0Nu || fCurrDecayMode == kHNLDcyNuNuNu
301 || fCurrDecayMode == kHNLDcyPi0Pi0Nu || fCurrDecayMode == kHNLDcyTEST );
302 }
303 //assert( event->FinalStatePrimaryLepton() );
304
305 // loop over all FS particles and set their mother to HNL
306 int itmp = 1, ilast = 1;
307 while( event->Particle( itmp ) ){
308 if( event->Particle(itmp)->Status() != kIStStableFinalState ){ itmp++; continue; }
309 event->Particle(itmp)->SetFirstMother(0);
310 event->Particle(itmp)->SetLastMother(-1);
311 if( itmp != ilead ) ilast = itmp;
312 itmp++;
313 }
314 event->Particle( 0 )->SetLastDaughter( ilast );
315 // "last daughter" of HNL means last non-primary-FS-lepton, so can be less than "first" daughter
316
317 LOG("HNL", pNOTICE)
318 << "Finished with decay products. Clean up and exit!";
319
320 // Clean-up
321 delete [] mass;
322 delete p4d;
323 delete v4d;
324}
pNOTICE
#define pNOTICE
Definition Messenger.h:61
pINFO
#define pINFO
Definition Messenger.h:62
pERROR
#define pERROR
Definition Messenger.h:59
genie::GHepParticle::GetP4
TLorentzVector * GetP4(void) const
Definition GHepParticle.cxx:194
genie::GHepParticle::GetX4
TLorentzVector * GetX4(void) const
Definition GHepParticle.cxx:212
genie::GHepParticle::E
double E(void) const
Get energy.
Definition GHepParticle.h:91
genie::GHepParticle::Status
GHepStatus_t Status(void) const
Definition GHepParticle.h:64
genie::GHepRecord::Probe
virtual GHepParticle * Probe(void) const
Definition GHepRecord.cxx:284
genie::GHepRecord::FinalStatePrimaryLepton
virtual GHepParticle * FinalStatePrimaryLepton(void) const
Definition GHepRecord.cxx:333
genie::Interaction::ExclTagPtr
XclsTag * ExclTagPtr(void) const
Definition Interaction.h:77
genie::XclsTag::SetNPions
void SetNPions(int npi_plus, int npi_0, int npi_minus)
Definition XclsTag.cxx:88
genie::hnl::Decayer::PolarisedDecay
bool PolarisedDecay(TGenPhaseSpace &fPSG, PDGCodeList pdgv, double wm, TVector3 vPolDir) const
Definition HNLDecayer.cxx:590
genie::hnl::Decayer::fCurrDecayMode
genie::hnl::HNLDecayMode_t fCurrDecayMode
Definition HNLDecayer.h:100
genie::hnl::Decayer::fPolDir
std::vector< double > fPolDir
Definition HNLDecayer.h:103
genie::hnl::Decayer::UnpolarisedDecay
bool UnpolarisedDecay(TGenPhaseSpace &fPSG, PDGCodeList pdgv, double wm) const
Definition HNLDecayer.cxx:550
genie::utils::hnl::DecayProductList
PDGCodeList DecayProductList(genie::hnl::HNLDecayMode_t hnldm)
genie::utils::print::P4AsString
string P4AsString(const TLorentzVector *p)
Definition PrintUtils.cxx:27
genie::kPdgPiM
const int kPdgPiM
Definition PDGCodes.h:159
genie::kIStStableFinalState
@ kIStStableFinalState
Definition GHepStatus.h:30
genie::kPdgPi0
const int kPdgPi0
Definition PDGCodes.h:160
genie::GHepStatus_t
enum genie::EGHepStatus GHepStatus_t
genie::kPdgPiP
const int kPdgPiP
Definition PDGCodes.h:158

References genie::utils::hnl::DecayProductList(), genie::GHepParticle::E(), genie::Interaction::ExclTagPtr(), fCurrDecayMode, fCurrInitStatePdg, fDecHadPdg, fDecLepPdg, fDoPol, fGetCMFrameInstead, genie::GHepRecord::FinalStatePrimaryLepton(), genie::PDGLibrary::Find(), fIsMajorana, fPolDir, genie::PDGLibrary::Instance(), genie::RandomGen::Instance(), genie::hnl::kHNLDcyNuNuNu, genie::hnl::kHNLDcyPi0Nu, genie::hnl::kHNLDcyPi0Pi0Nu, genie::hnl::kHNLDcyTEST, genie::kIStStableFinalState, genie::kPdgPi0, genie::kPdgPiM, genie::kPdgPiP, LOG, genie::utils::print::P4AsString(), genie::GHepRecord::Particle(), pDEBUG, pERROR, pINFO, pNOTICE, PolarisedDecay(), genie::GHepRecord::Probe(), genie::PDGCodeList::push_back(), pWARN, genie::RandomGen::RndGen(), genie::XclsTag::SetNPions(), genie::exceptions::EVGThreadException::SetReason(), genie::GHepParticle::Status(), genie::exceptions::EVGThreadException::SwitchOnFastForward(), and UnpolarisedDecay().

Referenced by ProcessEventRecord().

◆ GenerateMomentum()

std::vector< double > * Decayer::GenerateMomentum ( GHepRecord * event) const
private

Definition at line 351 of file HNLDecayer.cxx.

352{
353 Interaction * interaction = event->Summary();
354 double E = interaction->InitState().ProbeE(kRfLab);
355 double M = PDGLibrary::Instance()->Find(kPdgHNL)->Mass();
356 double p = TMath::Sqrt(E*E-M*M);
357
358 // set some initial deviation from beam axis due to collimation effect
359 double thetaDev = 0.0; //fAngularDeviation; // deg
360 thetaDev *= genie::constants::kPi / 180.0; // rad
361 RandomGen * Rng = RandomGen::Instance();
362 double theta = Rng->RndGen().Gaus(0.0, thetaDev);
363 if( theta < 0.0 ) theta *= -1.0;
364 double phi = Rng->RndGen().Uniform(0.0, 2.0 * genie::constants::kPi);
365
366 double px = p * std::sin(theta) * std::cos(phi);
367 double py = p * std::sin(theta) * std::sin(phi);
368 double pz = p * std::cos(theta);
369
370 std::vector< double > * p3HNL = new std::vector< double >();
371 p3HNL->emplace_back(px);
372 p3HNL->emplace_back(py);
373 p3HNL->emplace_back(pz);
374
375 return p3HNL;
376}
genie::kPdgHNL
const int kPdgHNL
Definition PDGCodes.h:224

References genie::PDGLibrary::Find(), genie::Interaction::InitState(), genie::PDGLibrary::Instance(), genie::RandomGen::Instance(), genie::kPdgHNL, genie::constants::kPi, genie::kRfLab, genie::InitialState::ProbeE(), and genie::RandomGen::RndGen().

Referenced by AddInitialState().

◆ GetHNLCouplings()

std::vector< double > Decayer::GetHNLCouplings ( ) const
virtual

Implements genie::hnl::DecayRecordVisitorI.

Definition at line 715 of file HNLDecayer.cxx.

716{
717 std::vector<double> allCoups;
718 allCoups.emplace_back(fUe42); allCoups.emplace_back(fUm42); allCoups.emplace_back(fUt42);
719 return allCoups;
720}

References fUe42, fUm42, and fUt42.

Referenced by main(), and ReadInConfig().

◆ GetHNLInstance()

SimpleHNL Decayer::GetHNLInstance ( ) const
private

Definition at line 499 of file HNLDecayer.cxx.

500{
501 SimpleHNL sh = SimpleHNL( "HNLInstance", 0, genie::kPdgHNL, genie::kPdgKP,
502 fMass, fUe42, fUm42, fUt42, fIsMajorana );
503 sh.SetType( 0 );
504 sh.SetInterestingChannelsVec( fIntChannels );
505 sh.SetAngularDeviation( 0.0 ); //fAngularDeviation );
506 sh.SetBeam2UserTranslation( fTx, fTy, fTz );
507 sh.SetBeam2UserRotation( fR1, fR2, fR3 );
508 return sh;
509}
genie::hnl::Decayer::fIntChannels
std::vector< genie::hnl::HNLDecayMode_t > fIntChannels
Definition HNLDecayer.h:114
genie::hnl::SimpleHNL::SetBeam2UserTranslation
void SetBeam2UserTranslation(const double tx, const double ty, const double tz)
Definition SimpleHNL.h:298
genie::hnl::SimpleHNL::SetBeam2UserRotation
void SetBeam2UserRotation(const double r1, const double r2, const double r3)
Definition SimpleHNL.h:302
genie::hnl::SimpleHNL::SetInterestingChannelsVec
void SetInterestingChannelsVec(const std::vector< genie::hnl::HNLDecayMode_t > decVec)
Definition SimpleHNL.h:281
genie::hnl::SimpleHNL::SetAngularDeviation
void SetAngularDeviation(const double adev)
Definition SimpleHNL.h:296
genie::hnl::SimpleHNL::SetType
void SetType(const int type)
Definition SimpleHNL.h:294
genie::kPdgKP
const int kPdgKP
Definition PDGCodes.h:172

References fIntChannels, fIsMajorana, fMass, fR1, fR2, fR3, fTx, fTy, fTz, fUe42, fUm42, fUt42, genie::kPdgHNL, genie::kPdgKP, genie::hnl::SimpleHNL::SetAngularDeviation(), genie::hnl::SimpleHNL::SetBeam2UserRotation(), genie::hnl::SimpleHNL::SetBeam2UserTranslation(), genie::hnl::SimpleHNL::SetInterestingChannelsVec(), and genie::hnl::SimpleHNL::SetType().

Referenced by GetHNLLifetime(), and LoadConfig().

◆ GetHNLInterestingChannels()

std::string Decayer::GetHNLInterestingChannels ( ) const
virtual

Implements genie::hnl::DecayRecordVisitorI.

Definition at line 727 of file HNLDecayer.cxx.

728{
729 std::string chanInt = "";
730 for( int iCBits = sizeof( fChanBits ) / sizeof( fChanBits[0] ) - 1; iCBits >= 0 ; iCBits-- ){
731 chanInt.append( Form("%d", fChanBits[iCBits]) );
732 }
733 return chanInt;
734}

References fChanBits.

Referenced by main(), and ReadInConfig().

◆ GetHNLLifetime()

double Decayer::GetHNLLifetime ( ) const
virtual

Implements genie::hnl::DecayRecordVisitorI.

Definition at line 705 of file HNLDecayer.cxx.

706{
707 return (this->GetHNLInstance()).GetCoMLifetime();
708}
genie::hnl::Decayer::GetHNLInstance
genie::hnl::SimpleHNL GetHNLInstance() const
Definition HNLDecayer.cxx:499

References GetHNLInstance().

Referenced by main(), and ReadInConfig().

◆ GetHNLMass()

double Decayer::GetHNLMass ( ) const
virtual

Implements genie::hnl::DecayRecordVisitorI.

Definition at line 710 of file HNLDecayer.cxx.

711{
712 return fMass;
713}

References fMass.

Referenced by GetCommandLineArgs(), main(), and ReadInConfig().

◆ GetPGunDeviation()

std::vector< double > Decayer::GetPGunDeviation ( ) const
virtual

Implements genie::hnl::DecayRecordVisitorI.

Definition at line 762 of file HNLDecayer.cxx.

763{
764 std::vector< double > PGDeviation;
765 PGDeviation.emplace_back( fPGDTheta ); PGDeviation.emplace_back( fPGDPhi );
766 return PGDeviation;
767}

References fPGDPhi, and fPGDTheta.

Referenced by GenerateOriginMomentum().

◆ GetPGunDirection()

std::vector< double > Decayer::GetPGunDirection ( ) const
virtual

Implements genie::hnl::DecayRecordVisitorI.

Definition at line 755 of file HNLDecayer.cxx.

756{
757 std::vector< double > PGDirection;
758 PGDirection.emplace_back( fPGCx ); PGDirection.emplace_back( fPGCy ); PGDirection.emplace_back( fPGCz );
759 return PGDirection;
760}

References fPGCx, fPGCy, and fPGCz.

Referenced by GenerateOriginMomentum(), and ReadInConfig().

◆ GetPGunDOrigin()

std::vector< double > Decayer::GetPGunDOrigin ( ) const
virtual

Implements genie::hnl::DecayRecordVisitorI.

Definition at line 743 of file HNLDecayer.cxx.

744{
745 std::vector< double > PGDOrigin;
746 PGDOrigin.emplace_back( fPGDx ); PGDOrigin.emplace_back( fPGDy ); PGDOrigin.emplace_back( fPGDz );
747 return PGDOrigin;
748}

References fPGDx, fPGDy, and fPGDz.

Referenced by GenerateOriginPosition().

◆ GetPGunEnergy()

double Decayer::GetPGunEnergy ( ) const
virtual

Implements genie::hnl::DecayRecordVisitorI.

Definition at line 750 of file HNLDecayer.cxx.

751{
752 return fPGE;
753}

References fPGE.

Referenced by main(), ReadInConfig(), and TestDecay().

◆ GetPGunOrigin()

std::vector< double > Decayer::GetPGunOrigin ( ) const
virtual

Implements genie::hnl::DecayRecordVisitorI.

Definition at line 736 of file HNLDecayer.cxx.

737{
738 std::vector< double > PGOrigin;
739 PGOrigin.emplace_back( fPGOx ); PGOrigin.emplace_back( fPGOy ); PGOrigin.emplace_back( fPGOz );
740 return PGOrigin;
741}

References fPGOx, fPGOy, and fPGOz.

Referenced by GenerateOriginPosition(), and TestDecay().

◆ IsHNLMajorana()

bool Decayer::IsHNLMajorana ( ) const
virtual

Implements genie::hnl::DecayRecordVisitorI.

Definition at line 722 of file HNLDecayer.cxx.

723{
724 return fIsMajorana;
725}

References fIsMajorana.

Referenced by main(), and ReadInConfig().

◆ LoadConfig()

void Decayer::LoadConfig ( void )
private

Definition at line 419 of file HNLDecayer.cxx.

420{
421 LOG("HNL", pDEBUG)
422 << "Loading configuration from file...";
423
424 this->GetParam( "HNL-Mass", fMass );
425 std::vector< double > U4l2s;
426 this->GetParamVect( "HNL-LeptonMixing", U4l2s );
427 SetHNLCouplings( U4l2s.at(0), U4l2s.at(1), U4l2s.at(2) );
428 this->GetParam( "HNL-Majorana", fIsMajorana );
429 //this->GetParam( "HNL-Type", fType );
430
431 //this->GetParam( "HNL-angular_deviation", fAngularDeviation );
432
433 this->GetParam( "IncludePolarisation", fDoPol );
434
435 this->GetParamVect( "Near2User_T", fB2UTranslation );
436 this->GetParamVect( "Near2Beam_R", fB2URotation );
437 SetBeam2User( fB2UTranslation, fB2URotation );
438
439 fIntChannels = {}; bool itChan = false;
440 //int chanBits[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
441
442 this->GetParam( "HNL-3B_nu_nu_nu", itChan ); if( itChan ){ fIntChannels.push_back( kHNLDcyNuNuNu ); fChanBits[0] = 1; }
443 this->GetParam( "HNL-3B_nu_e_e", itChan ); if( itChan ){ fIntChannels.push_back( kHNLDcyNuEE ); fChanBits[1] = 1; }
444 this->GetParam( "HNL-3B_nu_mu_e", itChan ); if( itChan ){ fIntChannels.push_back( kHNLDcyNuMuE ); fChanBits[2] = 1; }
445 this->GetParam( "HNL-2B_nu_pi0", itChan ); if( itChan ){ fIntChannels.push_back( kHNLDcyPi0Nu ); fChanBits[3] = 1; }
446 this->GetParam( "HNL-2B_e_pi", itChan ); if( itChan ){ fIntChannels.push_back( kHNLDcyPiE ); fChanBits[4] = 1; }
447 this->GetParam( "HNL-3B_nu_mu_mu", itChan ); if( itChan ){ fIntChannels.push_back( kHNLDcyNuMuMu ); fChanBits[5] = 1; }
448 this->GetParam( "HNL-2B_mu_pi", itChan ); if( itChan ){ fIntChannels.push_back( kHNLDcyPiMu ); fChanBits[6] = 1; }
449 this->GetParam( "HNL-3B_nu_pi0_pi0", itChan ); if( itChan ){ fIntChannels.push_back( kHNLDcyPi0Pi0Nu ); fChanBits[7] = 1; }
450 this->GetParam( "HNL-3B_e_pi_pi0", itChan ); if( itChan ){ fIntChannels.push_back( kHNLDcyPiPi0E ); fChanBits[8] = 1; }
451 this->GetParam( "HNL-3B_mu_pi_pi0", itChan ); if( itChan ){ fIntChannels.push_back( kHNLDcyPiPi0Mu ); fChanBits[9] = 1; }
452
453 this->GetParam( "GetCMFrameInstead", fGetCMFrameInstead );
454
455 // call GetHNLInstance here, to get lifetime
456 SimpleHNL sh = this->GetHNLInstance();
457 double CoMLifetime = sh.GetCoMLifetime();
458 assert( CoMLifetime > 0.0 );
459
460 // also read in particle gun parameters
461 this->GetParam( "PG-OriginX", fPGOx );
462 this->GetParam( "PG-OriginY", fPGOy );
463 this->GetParam( "PG-OriginZ", fPGOz );
464
465 this->GetParam( "PG-OriginDX", fPGDx );
466 this->GetParam( "PG-OriginDY", fPGDy );
467 this->GetParam( "PG-OriginDZ", fPGDz );
468
469 this->GetParam( "PG-Energy", fPGE );
470
471 this->GetParam( "PG-cx", fPGCx );
472 this->GetParam( "PG-cy", fPGCy );
473 this->GetParam( "PG-cz", fPGCz );
474
475 this->GetParam( "PG-DTheta", fPGDTheta );
476 this->GetParam( "PG-DPhi", fPGDPhi );
477
478 fIsConfigLoaded = true;
479}
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
genie::Algorithm::GetParamVect
int GetParamVect(const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
Handle to load vectors of parameters.
genie::hnl::Decayer::SetHNLCouplings
void SetHNLCouplings(double Ue42, double Um42, double Ut42) const
Definition HNLDecayer.cxx:481
genie::hnl::Decayer::fB2URotation
std::vector< double > fB2URotation
Definition HNLDecayer.h:120
genie::hnl::Decayer::fB2UTranslation
std::vector< double > fB2UTranslation
Definition HNLDecayer.h:118
genie::hnl::Decayer::SetBeam2User
void SetBeam2User(std::vector< double > translation, std::vector< double > rotation) const
Definition HNLDecayer.cxx:488
genie::hnl::SimpleHNL::GetCoMLifetime
double GetCoMLifetime()
Definition SimpleHNL.h:96
CoMLifetime
double CoMLifetime
Definition gBeamHNLEvGen.cxx:212

References CoMLifetime, fB2URotation, fB2UTranslation, fChanBits, fDoPol, fGetCMFrameInstead, fIntChannels, fIsConfigLoaded, fIsMajorana, fMass, fPGCx, fPGCy, fPGCz, fPGDPhi, fPGDTheta, fPGDx, fPGDy, fPGDz, fPGE, fPGOx, fPGOy, fPGOz, genie::hnl::SimpleHNL::GetCoMLifetime(), GetHNLInstance(), genie::Algorithm::GetParam(), genie::Algorithm::GetParamVect(), genie::hnl::kHNLDcyNuEE, genie::hnl::kHNLDcyNuMuE, genie::hnl::kHNLDcyNuMuMu, genie::hnl::kHNLDcyNuNuNu, genie::hnl::kHNLDcyPi0Nu, genie::hnl::kHNLDcyPi0Pi0Nu, genie::hnl::kHNLDcyPiE, genie::hnl::kHNLDcyPiMu, genie::hnl::kHNLDcyPiPi0E, genie::hnl::kHNLDcyPiPi0Mu, LOG, pDEBUG, SetBeam2User(), and SetHNLCouplings().

Referenced by Configure(), and Configure().

◆ PolarisedDecay()

bool Decayer::PolarisedDecay ( TGenPhaseSpace & fPSG,
PDGCodeList pdgv,
double wm,
TVector3 vPolDir ) const
private

Definition at line 590 of file HNLDecayer.cxx.

591{
592 // calculate polarisation modulus
593 PDGLibrary * pdgl = PDGLibrary::Instance();
594 double MHNL = pdgl->Find( kPdgHNL )->Mass();
595 double polMag = this->CalcPolMag( fParentPdg, fProdLepPdg, MHNL );
596 double polMod = -999.99;
597
598 // do decays until weight \in Uniform[0.0, 2.0] < 1 \mp polMod * cos\theta
599 unsigned int iUPD = 0;
600 double polWgt = -999.9;
601 RandomGen * rnd = RandomGen::Instance();
602 double rwgt = 0.0;
603 bool isAccepted = false;
604
605 bool failed = false;
606
607 // first, check to see if we have pi0 + v. Then let the neutrino be a QLep.
608 bool isPi0Nu = ( pdgv.size() == 3 &&
609 std::abs( pdgv.at(1) ) == kPdgPi0 &&
610 std::abs( pdgv.at(2) ) == kPdgNuMu );
611
612 while( !isAccepted && iUPD < controls::kMaxUnweightDecayIterations ){
613 bool failed = this->UnpolarisedDecay( fPSG, pdgv, wm );
614
615 // find charged lepton of FS. If two, take the leading one.
616 // For now, this method doesn't handle vvv invisible decay mode.
617
618 TVector3 lepDir;
619 vector<int>::const_iterator pdg_iter;
620 int idc = 0; double Elead = -1.0;
621 for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
622 int pdgc = *pdg_iter; Elead = -1.0;
623 bool isQLep = ( std::abs( pdgc ) == kPdgElectron ||
624 std::abs( pdgc ) == kPdgMuon ||
625 std::abs( pdgc ) == kPdgTau ||
626 ( isPi0Nu && std::abs( pdgc ) == kPdgNuMu ) );
627 if( isQLep ){
628 TLorentzVector * p4lep = fPSG.GetDecay(idc);
629 if( p4lep->E() > Elead ){
630 Elead = p4lep->E();
631 lepDir = p4lep->Vect();
632 fDecLepPdg = pdgc;
633 if( std::abs(fDecLepPdg) != std::abs(pdgc) || polMod < -1.0 ){
634 // update polarisation modulus for new leading lepton
635 polMod = this->CalcPolMod( polMag, fDecLepPdg, fDecHadPdg, MHNL );
636 } // std::abs(fDecLepPdg) != std::abs(pdgc) || polMod == -999.9
637 } // p4lep->E() > Elead
638 } // isQLep
639
640 idc++;
641 }
642
643 // find angle \theta of leading FS QLep with vPolDir
644 // assume differential decay rate \propto ( 1 \pm pMod * cos\theta )
645
646 double theta = vPolDir.Angle( lepDir ); // rad
647 double ctheta = TMath::Cos( theta );
648
649 rwgt = rnd->RndGen().Uniform(0.0, 2.0);
650 int typeMod = ( *(pdgv.begin()) > 0 ) ? 1 : -1;
651 polWgt = 1 + typeMod * polMod * ctheta;
652
653 isAccepted = ( rwgt >= polWgt );
654
655 iUPD++;
656 } // while( rwgt >= polWgt && iUPD < controls::kMaxUnweightDecayIterations )
657
658 if( iUPD == controls::kMaxUnweightDecayIterations ){
659 // report and return
660 LOG("HNL", pWARN)
661 << "Couldn't generate a polarised decay after "
662 << iUPD << " attempts";
663 failed = true;
664 return failed;
665 }
666
667 return failed;
668
669}
genie::hnl::Decayer::CalcPolMag
double CalcPolMag(int parPdg, int lepPdg, double M) const
Definition HNLDecayer.cxx:671
genie::hnl::Decayer::CalcPolMod
double CalcPolMod(double polMag, int lepPdg, int hadPdg, double M) const
Definition HNLDecayer.cxx:688
genie::controls::kMaxUnweightDecayIterations
static const unsigned int kMaxUnweightDecayIterations
Definition Controls.h:61
genie::kPdgTau
const int kPdgTau
Definition PDGCodes.h:39
genie::kPdgMuon
const int kPdgMuon
Definition PDGCodes.h:37
genie::kPdgNuMu
const int kPdgNuMu
Definition PDGCodes.h:30
genie::kPdgElectron
const int kPdgElectron
Definition PDGCodes.h:35

References CalcPolMag(), CalcPolMod(), fDecHadPdg, fDecLepPdg, genie::PDGLibrary::Find(), fParentPdg, fProdLepPdg, genie::PDGLibrary::Instance(), genie::RandomGen::Instance(), genie::controls::kMaxUnweightDecayIterations, genie::kPdgElectron, genie::kPdgHNL, genie::kPdgMuon, genie::kPdgNuMu, genie::kPdgPi0, genie::kPdgTau, LOG, pWARN, genie::RandomGen::RndGen(), and UnpolarisedDecay().

Referenced by GenerateDecayProducts().

◆ ProcessEventRecord()

void Decayer::ProcessEventRecord ( GHepRecord * event) const
virtual

Implements genie::hnl::DecayRecordVisitorI.

Definition at line 54 of file HNLDecayer.cxx.

55{
56
57 Interaction * interaction = event->Summary();
58
59 fCurrInitStatePdg = interaction->InitState().ProbePdg();
60 fCurrDecayMode = (HNLDecayMode_t) interaction->ExclTag().DecayMode();
61
62 LOG("HNL", pNOTICE)
63 << "Simulating HNL decay " << utils::hnl::AsString(fCurrDecayMode)
64 << " for an initial state with PDG code " << fCurrInitStatePdg;
65
66 this->ReadCreationInfo(event);
67 this->AddInitialState(event);
68 this->GenerateDecayProducts(event);
69 this->UpdateEventRecord(event);
70
71}
genie::Interaction::ExclTag
const XclsTag & ExclTag(void) const
Definition Interaction.h:72
genie::XclsTag::DecayMode
int DecayMode(void) const
Definition XclsTag.h:70
genie::hnl::Decayer::UpdateEventRecord
void UpdateEventRecord(GHepRecord *event) const
Definition HNLDecayer.cxx:378
genie::hnl::Decayer::AddInitialState
void AddInitialState(GHepRecord *event) const
Definition HNLDecayer.cxx:73
genie::hnl::Decayer::GenerateDecayProducts
void GenerateDecayProducts(GHepRecord *event) const
Definition HNLDecayer.cxx:121
genie::hnl::Decayer::ReadCreationInfo
void ReadCreationInfo(GHepRecord *event) const
Definition HNLDecayer.cxx:518
genie::hnl::HNLDecayMode_t
enum genie::hnl::EHNLDecayMode HNLDecayMode_t
genie::utils::hnl::AsString
string AsString(genie::hnl::HNLDecayMode_t hnldm)

References AddInitialState(), genie::utils::hnl::AsString(), genie::XclsTag::DecayMode(), genie::Interaction::ExclTag(), fCurrDecayMode, fCurrInitStatePdg, GenerateDecayProducts(), genie::Interaction::InitState(), LOG, pNOTICE, genie::InitialState::ProbePdg(), ReadCreationInfo(), and UpdateEventRecord().

Referenced by main(), and TestDecay().

◆ ReadCreationInfo()

void Decayer::ReadCreationInfo ( GHepRecord * event) const
private

Definition at line 518 of file HNLDecayer.cxx.

519{
520 if( fPolDir.size() > 0 ) fPolDir.clear();
521 /*
522 fPolDir.emplace_back( gnmf.ppvx );
523 fPolDir.emplace_back( gnmf.ppvy );
524 fPolDir.emplace_back( gnmf.ppvz );
525
526 fParentPdg = gnmf.ptype;
527 fProdLepPdg = gnmf.ppmedium;
528 */
529
530 TLorentzVector * vv = event->Vertex();
531 TLorentzVector * tmpx4 = event->Particle(0)->GetX4();
532 // this will only have been modified if we have enough polarisation information to do something.
533 // Otherwise, FluxCreator hasn't been run and we shouldn't do stuff.
534 if( tmpx4->X() != vv->X() || tmpx4->Y() != vv->Y() || tmpx4->Z() != vv->Z() ){
535 fPolDir.emplace_back( tmpx4->X() );
536 fPolDir.emplace_back( tmpx4->Y() );
537 int tmpMod = ( event->XSec() > 0.0 ) ? 1 : -1;
538 fPolDir.emplace_back( tmpMod * std::sqrt( 1.0 -
539 ( tmpx4->X() * tmpx4->X() + tmpx4->Y() * tmpx4->Y() ) ) );
540
541 fParentPdg = tmpx4->Z();
542 fProdLepPdg = tmpx4->T();
543
544 // now reset the x4 of the HNL to whatever the vertex is
545 event->Particle(0)->SetPosition( vv->X(), vv->Y(), vv->Z(), vv->T() );
546 event->SetXSec(0.0);
547 } else { return; }
548}

References fParentPdg, fPolDir, and fProdLepPdg.

Referenced by ProcessEventRecord().

◆ SetBeam2User()

void Decayer::SetBeam2User ( std::vector< double > translation,
std::vector< double > rotation ) const
private

Definition at line 488 of file HNLDecayer.cxx.

489{
490 fTx = -1.0 * translation.at(0);
491 fTy = -1.0 * translation.at(1);
492 fTz = -1.0 * translation.at(2);
493
494 fR1 = rotation.at(0);
495 fR2 = rotation.at(1);
496 fR3 = rotation.at(2);
497}

References fR1, fR2, fR3, fTx, fTy, and fTz.

Referenced by LoadConfig().

◆ SetHNLCouplings()

void Decayer::SetHNLCouplings ( double Ue42,
double Um42,
double Ut42 ) const
private

Definition at line 481 of file HNLDecayer.cxx.

482{
483 fUe42 = Ue42;
484 fUm42 = Um42;
485 fUt42 = Ut42;
486}

References fUe42, fUm42, and fUt42.

Referenced by LoadConfig().

◆ SetProdVtxPosition()

void Decayer::SetProdVtxPosition ( const TLorentzVector & v4) const
private

Definition at line 511 of file HNLDecayer.cxx.

512{
513 TLorentzVector * pv4 = new TLorentzVector();
514 pv4->SetXYZT( v4.X(), v4.Y(), v4.Z(), v4.T() );
515 fProdVtx = pv4;
516}
genie::hnl::Decayer::fProdVtx
TLorentzVector * fProdVtx
Definition HNLDecayer.h:124

References fProdVtx.

Referenced by GenerateDecayPosition().

◆ UnpolarisedDecay()

bool Decayer::UnpolarisedDecay ( TGenPhaseSpace & fPSG,
PDGCodeList pdgv,
double wm ) const
private

accept_decay

Definition at line 550 of file HNLDecayer.cxx.

551{
552
553 RandomGen * rnd = RandomGen::Instance();
554
555 bool accept_decay=false;
556 unsigned int itry=0;
557
558 bool failed = false;
559 while(!accept_decay) {
560 itry++;
561
562 if(itry > controls::kMaxUnweightDecayIterations) {
563 // report and return
564 LOG("HNL", pWARN)
565 << "Couldn't generate an unweighted phase space decay after "
566 << itry << " attempts";
567 failed = true;
568 return failed;
569 }
570 double w = fPSG.Generate();
571 if(w > wm) {
572 LOG("HNL", pWARN)
573 << "Decay weight = " << w << " > max decay weight = " << wm;
574 }
575 double gw = wm * rnd->RndHadro().Rndm();
576 accept_decay = (gw<=w);
577
578 /*
579 LOG("HNL", pDEBUG)
580 << "Decay weight = " << w << " / R = " << gw
581 << " - accepted: " << accept_decay;
582 */
583
584 } //!accept_decay
585
586 return failed;
587
588}
genie::RandomGen::RndHadro
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition RandomGen.h:53

References genie::RandomGen::Instance(), genie::controls::kMaxUnweightDecayIterations, LOG, pWARN, and genie::RandomGen::RndHadro().

Referenced by GenerateDecayProducts(), and PolarisedDecay().

◆ UpdateEventRecord()

void Decayer::UpdateEventRecord ( GHepRecord * event) const
private

Definition at line 378 of file HNLDecayer.cxx.

379{
380 Interaction * interaction = event->Summary();
381
382 interaction->KinePtr()->Sett( 0.0, true );
383 interaction->KinePtr()->SetW( interaction->InitStatePtr()->Probe()->Mass(), true );
384 TLorentzVector * p4HNL = event->Particle(0)->GetP4(); assert( p4HNL );
385 // primary lepton is FirstDaughter() of Probe()
386 // need Probe() as a GHepParticle(), not a TParticlePDG()!
387 // get from event record position 0
388 TLorentzVector * p4FSL = 0;
389 if( event->FinalStatePrimaryLepton() ){
390 int iFSL = event->Particle(0)->FirstDaughter();
391 assert( event->Particle( iFSL ) );
392 p4FSL = ( event->Particle( iFSL ) )->GetP4();
393 assert( p4FSL );
394 TLorentzVector p4DIF( p4HNL->Px() - p4FSL->Px(),
395 p4HNL->Py() - p4FSL->Py(),
396 p4HNL->Pz() - p4FSL->Pz(),
397 p4HNL->E() - p4FSL->E() );
398 interaction->KinePtr()->SetQ2( p4DIF.M2(), true );
399
400 }
401
402 // clean up
403 delete p4HNL;
404 if(p4FSL) delete p4FSL;
405}
genie::InitialState::Probe
TParticlePDG * Probe(void) const
Definition InitialState.cxx:147
genie::Interaction::KinePtr
Kinematics * KinePtr(void) const
Definition Interaction.h:76
genie::Kinematics::SetQ2
void SetQ2(double Q2, bool selected=false)
Definition Kinematics.cxx:255
genie::Kinematics::Sett
void Sett(double t, bool selected=false)
Definition Kinematics.cxx:291
genie::Kinematics::SetW
void SetW(double W, bool selected=false)
Definition Kinematics.cxx:279

References genie::GHepRecord::FinalStatePrimaryLepton(), genie::Interaction::InitStatePtr(), genie::Interaction::KinePtr(), genie::GHepRecord::Particle(), genie::InitialState::Probe(), genie::Kinematics::SetQ2(), genie::Kinematics::Sett(), and genie::Kinematics::SetW().

Referenced by ProcessEventRecord().

Member Data Documentation

◆ fB2URotation

std::vector< double > genie::hnl::Decayer::fB2URotation
mutableprivate

Definition at line 120 of file HNLDecayer.h.

Referenced by LoadConfig().

◆ fB2UTranslation

std::vector< double > genie::hnl::Decayer::fB2UTranslation
mutableprivate

Definition at line 118 of file HNLDecayer.h.

Referenced by LoadConfig().

◆ fChanBits

int genie::hnl::Decayer::fChanBits[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
mutableprivate

Definition at line 115 of file HNLDecayer.h.

115{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };

Referenced by GetHNLInterestingChannels(), and LoadConfig().

◆ fCurrDecayMode

genie::hnl::HNLDecayMode_t genie::hnl::Decayer::fCurrDecayMode
mutableprivate

Definition at line 100 of file HNLDecayer.h.

Referenced by GenerateDecayProducts(), and ProcessEventRecord().

◆ fCurrInitStatePdg

int genie::hnl::Decayer::fCurrInitStatePdg
mutableprivate

Definition at line 99 of file HNLDecayer.h.

Referenced by GenerateDecayProducts(), and ProcessEventRecord().

◆ fDecHadPdg

int genie::hnl::Decayer::fDecHadPdg
private

Definition at line 102 of file HNLDecayer.h.

Referenced by GenerateDecayProducts(), and PolarisedDecay().

◆ fDecLepPdg

int genie::hnl::Decayer::fDecLepPdg
private

Definition at line 102 of file HNLDecayer.h.

Referenced by GenerateDecayProducts(), and PolarisedDecay().

◆ fDoPol

bool genie::hnl::Decayer::fDoPol = false
mutableprivate

Definition at line 112 of file HNLDecayer.h.

Referenced by GenerateDecayProducts(), and LoadConfig().

◆ fGetCMFrameInstead

bool genie::hnl::Decayer::fGetCMFrameInstead = false
mutableprivate

Definition at line 133 of file HNLDecayer.h.

Referenced by GenerateDecayProducts(), and LoadConfig().

◆ fIntChannels

std::vector< genie::hnl::HNLDecayMode_t > genie::hnl::Decayer::fIntChannels
mutableprivate

Definition at line 114 of file HNLDecayer.h.

Referenced by GetHNLInstance(), and LoadConfig().

◆ fIsConfigLoaded

bool genie::hnl::Decayer::fIsConfigLoaded = false
mutableprivate

Definition at line 105 of file HNLDecayer.h.

Referenced by LoadConfig().

◆ fIsMajorana

bool genie::hnl::Decayer::fIsMajorana = false
mutableprivate

Definition at line 109 of file HNLDecayer.h.

Referenced by GenerateDecayProducts(), GetHNLInstance(), IsHNLMajorana(), and LoadConfig().

◆ fISMom

TLorentzVector* genie::hnl::Decayer::fISMom = 0
mutableprivate

Definition at line 125 of file HNLDecayer.h.

◆ fMass

double genie::hnl::Decayer::fMass
mutableprivate

Definition at line 107 of file HNLDecayer.h.

Referenced by GetHNLInstance(), GetHNLMass(), and LoadConfig().

◆ fParentPdg

int genie::hnl::Decayer::fParentPdg
mutableprivate

Definition at line 102 of file HNLDecayer.h.

Referenced by PolarisedDecay(), and ReadCreationInfo().

◆ fPGCx

double genie::hnl::Decayer::fPGCx = 0.0
mutableprivate

Definition at line 130 of file HNLDecayer.h.

Referenced by GetPGunDirection(), and LoadConfig().

◆ fPGCy

double genie::hnl::Decayer::fPGCy = 0.0
private

Definition at line 130 of file HNLDecayer.h.

Referenced by GetPGunDirection(), and LoadConfig().

◆ fPGCz

double genie::hnl::Decayer::fPGCz = 0.0
private

Definition at line 130 of file HNLDecayer.h.

Referenced by GetPGunDirection(), and LoadConfig().

◆ fPGDPhi

double genie::hnl::Decayer::fPGDPhi = 0.0
private

Definition at line 131 of file HNLDecayer.h.

Referenced by GetPGunDeviation(), and LoadConfig().

◆ fPGDTheta

double genie::hnl::Decayer::fPGDTheta = 0.0
mutableprivate

Definition at line 131 of file HNLDecayer.h.

Referenced by GetPGunDeviation(), and LoadConfig().

◆ fPGDx

double genie::hnl::Decayer::fPGDx = 0.0
mutableprivate

Definition at line 128 of file HNLDecayer.h.

Referenced by GetPGunDOrigin(), and LoadConfig().

◆ fPGDy

double genie::hnl::Decayer::fPGDy = 0.0
private

Definition at line 128 of file HNLDecayer.h.

Referenced by GetPGunDOrigin(), and LoadConfig().

◆ fPGDz

double genie::hnl::Decayer::fPGDz = 0.0
private

Definition at line 128 of file HNLDecayer.h.

Referenced by GetPGunDOrigin(), and LoadConfig().

◆ fPGE

double genie::hnl::Decayer::fPGE = 0.0
mutableprivate

Definition at line 129 of file HNLDecayer.h.

Referenced by GetPGunEnergy(), and LoadConfig().

◆ fPGOx

double genie::hnl::Decayer::fPGOx = 0.0
mutableprivate

Definition at line 127 of file HNLDecayer.h.

Referenced by GetPGunOrigin(), and LoadConfig().

◆ fPGOy

double genie::hnl::Decayer::fPGOy = 0.0
private

Definition at line 127 of file HNLDecayer.h.

Referenced by GetPGunOrigin(), and LoadConfig().

◆ fPGOz

double genie::hnl::Decayer::fPGOz = 0.0
private

Definition at line 127 of file HNLDecayer.h.

Referenced by GetPGunOrigin(), and LoadConfig().

◆ fPolDir

std::vector<double> genie::hnl::Decayer::fPolDir
mutableprivate

Definition at line 103 of file HNLDecayer.h.

Referenced by GenerateDecayProducts(), and ReadCreationInfo().

◆ fProdLepPdg

int genie::hnl::Decayer::fProdLepPdg
private

Definition at line 102 of file HNLDecayer.h.

Referenced by PolarisedDecay(), and ReadCreationInfo().

◆ fProdVtx

TLorentzVector* genie::hnl::Decayer::fProdVtx = 0
mutableprivate

Definition at line 124 of file HNLDecayer.h.

Referenced by SetProdVtxPosition().

◆ fProdVtxHist

TH3D* genie::hnl::Decayer::fProdVtxHist = 0
mutableprivate

Definition at line 123 of file HNLDecayer.h.

Referenced by GenerateDecayPosition().

◆ fR1

double genie::hnl::Decayer::fR1 = -1.0
mutableprivate

Definition at line 121 of file HNLDecayer.h.

Referenced by GetHNLInstance(), and SetBeam2User().

◆ fR2

double genie::hnl::Decayer::fR2 = -1.0
private

Definition at line 121 of file HNLDecayer.h.

Referenced by GetHNLInstance(), and SetBeam2User().

◆ fR3

double genie::hnl::Decayer::fR3 = -1.0
private

Definition at line 121 of file HNLDecayer.h.

Referenced by GetHNLInstance(), and SetBeam2User().

◆ fTx

double genie::hnl::Decayer::fTx = -1.0
mutableprivate

Definition at line 119 of file HNLDecayer.h.

Referenced by GetHNLInstance(), and SetBeam2User().

◆ fTy

double genie::hnl::Decayer::fTy = -1.0
private

Definition at line 119 of file HNLDecayer.h.

Referenced by GetHNLInstance(), and SetBeam2User().

◆ fTz

double genie::hnl::Decayer::fTz = -1.0
private

Definition at line 119 of file HNLDecayer.h.

Referenced by GetHNLInstance(), and SetBeam2User().

◆ fUe42

double genie::hnl::Decayer::fUe42 = -1.0
mutableprivate

Definition at line 108 of file HNLDecayer.h.

Referenced by GetHNLCouplings(), GetHNLInstance(), and SetHNLCouplings().

◆ fUm42

double genie::hnl::Decayer::fUm42 = -1.0
private

Definition at line 108 of file HNLDecayer.h.

Referenced by GetHNLCouplings(), GetHNLInstance(), and SetHNLCouplings().

◆ fUt42

double genie::hnl::Decayer::fUt42 = -1.0
private

Definition at line 108 of file HNLDecayer.h.

Referenced by GetHNLCouplings(), GetHNLInstance(), and SetHNLCouplings().

The documentation for this class was generated from the following files:
Generated on for GENIEGenerator by doxygen 1.14.0