GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DMEOutgoingDarkGenerator.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 
7  Costas Andreopoulos <c.andreopoulos \at cern.ch>
8  University of Liverpool
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Feb 09, 2009 - CA
14  Moved into the DME package from its previous location (EVGModules package)
15  @ Feb 12, 2013 - CA (code from Rosen Matev)
16  Add mass_electron^2 term in kinematical calculation.
17 */
18 //____________________________________________________________________________
19 
20 #include <TMath.h>
21 #include <TVector3.h>
22 
32 
33 using namespace genie;
34 using namespace genie::constants;
35 
36 //___________________________________________________________________________
38 OutgoingDarkGenerator("genie::DMEOutgoingDarkGenerator")
39 {
40 
41 }
42 //___________________________________________________________________________
44 OutgoingDarkGenerator("genie::DMEOutgoingDarkGenerator", config)
45 {
46 
47 }
48 //___________________________________________________________________________
50 {
51 
52 }
53 //___________________________________________________________________________
55 {
56 // This method generates the final state primary lepton for DME events
57 
58  Interaction * interaction = evrec->Summary();
59  const InitialState & init_state = interaction->InitState();
60 
61  // Get selected kinematics
62  double y = interaction->Kine().y(true);
63  assert(y>0 && y<1);
64 
65  // Final state primary lepton PDG code
66  int pdgc = interaction->FSPrimLeptonPdg();
67  assert(pdgc!=0);
68 
69  // Compute the neutrino and muon energy
70  double Ev = init_state.ProbeE(kRfLab);
71  double El = y*Ev;
72 
73  LOG("LeptonicVertex", pINFO)
74  << "Ev = " << Ev << ", y = " << y << ", -> El = " << El;
75 
76  // Compute the momentum transfer and scattering angle
77  double Ev2 = TMath::Power(Ev,2);
78  double El2 = TMath::Power(El,2);
79  double me = kElectronMass;
80  double ml = interaction->FSPrimLepton()->Mass();
81  double ml2 = TMath::Power(ml,2);
82  double pl = TMath::Sqrt(El2-ml2);
83  double pv = TMath::Sqrt(Ev2-ml2);
84 
85  assert(El2>=ml2&&Ev2>=ml2);
86 
87  double Q2 = 2*(Ev-El)*me;
88  double costh = (El*Ev - ml2 -0.5*Q2)/pl/pv;
89  double sinth = TMath::Sqrt( TMath::Max(0., 1-TMath::Power(costh,2.)) );
90 
91  LOG("LeptonicVertex", pNOTICE)
92  << "Q2 = " << Q2 << ", cos(theta) = " << costh;
93 
94  //warn about overflow in costheta and ignore it if it is small or abort
95  if( TMath::Abs(costh)>1 ) {
96  LOG("LeptonicVertex", pWARN)
97  << "El = " << El << ", Ev = " << Ev << ", cos(theta) = " << costh;
98  //if(TMath::Abs(costh)-1<0.3) costh = 1.0; //why?
99  }
100  assert(TMath::Abs(costh)<=1);
101 
102  // Compute the p components along and perpendicular the v direction
103  double plp = pl * costh; // p(//)
104  double plt = pl * sinth; // p(-|)
105 
106  LOG("LeptonicVertex", pNOTICE)
107  << "fsl: E = " << El << ", |p//| = " << plp << "[pT] = " << plt;
108 
109  // Randomize transverse components
110  RandomGen * rnd = RandomGen::Instance();
111  double phi = 2*kPi * rnd->RndLep().Rndm();
112  double pltx = plt * TMath::Cos(phi);
113  double plty = plt * TMath::Sin(phi);
114 
115  // Take a unit vector along the neutrino direction
116  TVector3 unit_nudir = evrec->Probe()->P4()->Vect().Unit();
117 
118  // Rotate lepton momentum vector from the reference frame (x'y'z') where
119  // {z':(neutrino direction), z'x':(theta plane)} to the LAB
120  TVector3 p3l(pltx,plty,plp);
121  p3l.RotateUz(unit_nudir);
122 
123  // Lepton 4-momentum in the LAB
124  TLorentzVector p4l(p3l,El);
125 
126  // Create a GHepParticle and add it to the event record
127  this->AddToEventRecord(evrec, pdgc, p4l);
128 
129  // Set final state lepton polarization
130  this->SetPolarization(evrec);
131 }
132 //___________________________________________________________________________
TRandom3 & RndLep(void) const
rnd number generator used by final state primary lepton generators
Definition: RandomGen.h:62
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
static const double kElectronMass
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
Abstract class. Is used to pass common implementation to concrete implementations of the EventRecordV...
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:284
double y(bool selected=false) const
Definition: Kinematics.cxx:112
Summary information for an interaction.
Definition: Interaction.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
virtual void SetPolarization(GHepRecord *ev) const
const Kinematics & Kine(void) const
Definition: Interaction.h:71
virtual void AddToEventRecord(GHepRecord *ev, int pdgc, const TLorentzVector &p4) const
#define pINFO
Definition: Messenger.h:62
void ProcessEventRecord(GHepRecord *event_rec) const
#define pWARN
Definition: Messenger.h:60
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
const InitialState & InitState(void) const
Definition: Interaction.h:69
#define pNOTICE
Definition: Messenger.h:61
double ProbeE(RefFrame_t rf) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Initial State information.
Definition: InitialState.h:48