GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AMNuGammaGenerator.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  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 
23 
24 using namespace genie;
25 using namespace genie::constants;
26 
27 //___________________________________________________________________________
29 EventRecordVisitorI("genie::AMNuGammaGenerator")
30 {
31 
32 }
33 //___________________________________________________________________________
35 EventRecordVisitorI("genie::AMNuGammaGenerator", config)
36 {
37 
38 }
39 //___________________________________________________________________________
41 {
42 
43 }
44 //___________________________________________________________________________
46 {
47  this->AddPhoton(evrec);
48  this->AddFinalStateNeutrino(evrec);
49  this->AddRecoilNucleon(evrec);
50 //this->AddTargetRemnant(evrec);
51 }
52 //___________________________________________________________________________
54 {
55 // Adding the final state photon
56 //
57  LOG("AMNuGammaGenerator", pINFO) << "Adding final state photon";
58 
60 
61  // Get boost vector for transforms between LAB <-> NRF (nucleon rest frame)
62  GHepParticle * nuc = evrec->HitNucleon();
63  const TLorentzVector & p4nuc_lab = *(nuc->P4());
64  TVector3 beta = p4nuc_lab.BoostVector();
65 
66  // Get the neutrino 4-momentum at the LAB
67  GHepParticle * nu = evrec->Probe();
68  const TLorentzVector & p4v_lab = *(nu->P4());
69 
70  // Get the neutrino 4-momentum at the NRF
71  TLorentzVector p4v_nrf = p4v_lab;
72  p4v_nrf.Boost(-1.*beta);
73 
74  // Generate the photon cos(theta) with respect to the neutrino direction
75  // (at NRF) and a uniform azimuthal angle phi
76  double costheta_gamma = -1.0 + 2.0 * rnd->RndKine().Rndm();
77  double phi_gamma = 2.0 * kPi * rnd->RndKine().Rndm();
78 
79  // Generate the fraction of the neutrino energy taken by the photon
80  double efrac_gamma = rnd->RndKine().Rndm();
81 
82  // Calculate the photon energy at the NRF
83  double Ev_nrf = p4v_nrf.Energy();
84  double Egamma_nrf = Ev_nrf * efrac_gamma;
85 
86  // Calculate the photon momentum components at a rotated NRF (NRF')
87  // where z is along the neutrino direction
88  double sintheta_gamma = TMath::Sqrt(1-TMath::Power(costheta_gamma,2));
89  double pgamma_nrf_p = Egamma_nrf * costheta_gamma; // p(//)
90  double pgamma_nrf_t = Egamma_nrf * sintheta_gamma; // p(-|)
91  double px = pgamma_nrf_t * TMath::Sin(phi_gamma);
92  double py = pgamma_nrf_t * TMath::Cos(phi_gamma);
93  double pz = pgamma_nrf_p;
94 
95  // Take a unit vector along the neutrino direction @ the NRF
96  TVector3 unit_nudir = p4v_nrf.Vect().Unit();
97 
98  // Rotate the photon momentum vector from NRF' -> NRF
99  TVector3 p3gamma_nrf(px,py,pz);
100  p3gamma_nrf.RotateUz(unit_nudir);
101 
102  // Get the photon 4-momentum back at the LAB
103  TLorentzVector p4gamma_lab(p3gamma_nrf, Egamma_nrf);
104  p4gamma_lab.Boost(beta);
105 
106  // Add the photon at the event record
107  const TLorentzVector & vtx = *(nu->X4());
108  GHepParticle p(kPdgGamma,kIStStableFinalState,0,-1,-1,-1,p4gamma_lab,vtx);
109  evrec->AddParticle(p);
110 }
111 //___________________________________________________________________________
113 {
114 // Adding the final state neutrino
115 // Just use 4-momentum conservation (init_neutrino = photon + final_neutrino)
116 
117  LOG("AMNuGammaGenerator", pINFO) << "Adding final state neutrino";
118 
119  GHepParticle * nu = evrec->Probe(); // incoming v
120  GHepParticle * gamma = evrec->Particle(nu->FirstDaughter()); // gamma
121  assert(nu);
122  assert(gamma);
123 
124  const TLorentzVector & vtx = *(nu->X4()); // vtx
125 
126  const TLorentzVector & p4nu_lab = *(nu->P4());
127  const TLorentzVector & p4gamma_lab = *(gamma->P4());
128  TLorentzVector p4_lab = p4nu_lab - p4gamma_lab;
129 
130  GHepParticle p(nu->Pdg(), kIStStableFinalState, 0,-1,-1,-1, p4_lab, vtx);
131  evrec->AddParticle(p);
132 }
133 //___________________________________________________________________________
135 {
136 // Adding the recoil nucleon.
137 // Recoil is negligible at the H3 model. However, for nuclear targets, the
138 // hit nucleon was off-the-mass-shell (bound). Doing some magic here to bring
139 // it back on-the-mass-shell so that it can be INTRANUKE'ed and appear in
140 // the final state. The nucleon will keep its original (Fermi) 3-momentum but
141 // take some energy back from the remnant nucleus. No such tweaking takes
142 // place for free nucleon targets.
143 
144  LOG("AMNuGammaGenerator", pINFO) << "Adding recoil nucleon";
145 
146  // Get the interaction target
147  GHepParticle * tgt_nucleus = evrec->TargetNucleus();
148  bool is_nuclear_target = (tgt_nucleus!=0);
149 
150  // Get the hit nucleon
151  GHepParticle * hitnuc = evrec->HitNucleon();
152  assert(hitnuc);
153 
154  // Get the hit nucleon pdg code (= recoil nucleon pdg code)
155  int pdgc = hitnuc->Pdg();
156 
157  // Get the hit nucleon 4-momentum (LAB)
158  const TLorentzVector & p4n = *(hitnuc->P4());
159  TLorentzVector p4(p4n);
160 
161  // Tweak the 4-momentum to bring the recoil nucleon on-the-mass-shell
162  // (for nuclear targets only)
163  if (is_nuclear_target) {
164  double p = p4.Vect().Mag();
165  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
166  double E = TMath::Sqrt(m*m+p*p);
167  p4.SetE(E);
168  }
169 
170  // Get the vtx position
171  GHepParticle * neutrino = evrec->Probe();
172  const TLorentzVector & vtx = *(neutrino->X4());
173 
174  // Add the recoil nucleon at the event record
175  GHepStatus_t ist = (is_nuclear_target) ?
177  int mom = evrec->HitNucleonPosition();
178 
179  LOG("AMNuGammaGenerator", pINFO)
180  << "Adding recoil baryon [pdgc = " << pdgc << "]";
181 
182  GHepParticle p(pdgc, ist, mom,-1,-1,-1, p4, vtx);
183 // double w = hitnuc->RemovalEnergy();
184 /// p.SetRemovalEnergy(w);
185  evrec->AddParticle(p);
186 }
187 //___________________________________________________________________________
189 {
190 // Add the remnant nuclear target at the GHEP record
191 
192  LOG("AMNuGammaGenerator", pINFO) << "Adding final state nucleus";
193 
194  // Skip for non nuclear targets
195  GHepParticle * nucleus = evrec->TargetNucleus();
196  if (!nucleus) {
197  LOG("AMNuGammaGenerator", pDEBUG)
198  << "No nucleus in the initial state - no remnant nucleus to add in the f/s";
199  return;
200  }
201 
202  // Compute A,Z for final state nucleus & get its PDG code and its mass
203  //GHepParticle * nucleon = evrec->HitNucleon();
204 
205  GHepParticle * hit_nucleon = evrec->HitNucleon(); // hit
206  GHepParticle * rec_nucleon = evrec->Particle(hit_nucleon->FirstDaughter()); // recoil
207 
208  assert(rec_nucleon);
209  int npdgc = rec_nucleon->Pdg();
210  bool is_p = pdg::IsProton(npdgc);
211  int A = nucleus->A();
212  int Z = nucleus->Z();
213  if (is_p) Z--;
214  A--;
215  int ipdgc = pdg::IonPdgCode(A, Z);
216  TParticlePDG * remnant = PDGLibrary::Instance()->Find(ipdgc);
217  if(!remnant) {
218  LOG("AMNuGammaGenerator", pFATAL)
219  << "No particle with [A = " << A << ", Z = " << Z
220  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
221  assert(remnant);
222  }
223 
224  // Figure out the remnant 4-momentum
225  double Mf = remnant->Mass();
226  double Mf2 = TMath::Power(Mf,2);
227  double px = -1.* rec_nucleon->Px();
228  double py = -1.* rec_nucleon->Py();
229  double pz = -1.* rec_nucleon->Pz();
230  double E = TMath::Sqrt(Mf2 + rec_nucleon->P4()->Vect().Mag2());
231  E += (hit_nucleon->P4()->E() - rec_nucleon->P4()->E());
232 
233  //-- Add the nucleus to the event record
234  LOG("AMNuGammaGenerator", pINFO)
235  << "Adding nucleus [A = " << A << ", Z = " << Z
236  << ", pdgc = " << ipdgc << "]";
237  int mom = evrec->TargetNucleusPosition();
238  evrec->AddParticle(
239  ipdgc,kIStStableFinalState, mom,-1,-1,-1, px,py,pz,E, 0,0,0,0);
240 }
241 //___________________________________________________________________________
int Z(void) const
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
int FirstDaughter(void) const
Definition: GHepParticle.h:68
#define pFATAL
Definition: Messenger.h:56
virtual int HitNucleonPosition(void) const
Definition: GHepRecord.cxx:418
void AddTargetRemnant(GHepRecord *event_rec) const
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
void AddPhoton(GHepRecord *event_rec) const
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:284
int Pdg(void) const
Definition: GHepParticle.h:63
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double A
Definition: Units.h:74
const int kPdgGamma
Definition: PDGCodes.h:189
void AddFinalStateNeutrino(GHepRecord *event_rec) const
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
#define pINFO
Definition: Messenger.h:62
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:50
void AddRecoilNucleon(GHepRecord *event_rec) const
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:313
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
int A(void) const
void ProcessEventRecord(GHepRecord *event_rec) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
static constexpr double m
Definition: Units.h:71
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
virtual int TargetNucleusPosition(void) const
Definition: GHepRecord.cxx:370
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63
double Py(void) const
Get Py.
Definition: GHepParticle.h:89