AMNuGammaGenerator.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab - October 03, 2004
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Feb 25, 2008 - CA
14  This event generation modules was first added in version 2.3.1 as part of
15  the new event generation thread handling amonaly-mediated single gamma
16  interactions.
17  @ Sep 25, 2008 - CA
18  Improved the calculation of the photon 4-momentum. Generating the 4-momentum
19  at NRF' and then rotating it to NRF and boosting it back at the LAB.
20  The photon cos(theta) follows a uniform distribution (with respect to the
21  incoming neutrnino in NRF). Need further inputs for modeling the energy
22  transfer to the photon (uniform distribution up to the available energy).
23  The final state neutrino 4-momentum determined from energy conservation
24  (incoming nu = outgoing gamma + outgoing nu) but difficult to keep on the
25  mass shell. The nucleon recoil is negligible. For nuclear targets the hit
26  nucleon is forced back on the mass shell before intranuclear rescattering.
27  The necessary energy is taken from the remnant nucleus.
28 */
29 //____________________________________________________________________________
30 
31 #include <TMath.h>
32 
43 
44 using namespace genie;
45 using namespace genie::constants;
46 
47 //___________________________________________________________________________
49 EventRecordVisitorI("genie::AMNuGammaGenerator")
50 {
51 
52 }
53 //___________________________________________________________________________
55 EventRecordVisitorI("genie::AMNuGammaGenerator", config)
56 {
57 
58 }
59 //___________________________________________________________________________
61 {
62 
63 }
64 //___________________________________________________________________________
66 {
67  this->AddPhoton(evrec);
68  this->AddFinalStateNeutrino(evrec);
69  this->AddRecoilNucleon(evrec);
70 //this->AddTargetRemnant(evrec);
71 }
72 //___________________________________________________________________________
74 {
75 // Adding the final state photon
76 //
77  LOG("AMNuGammaGenerator", pINFO) << "Adding final state photon";
78 
80 
81  // Get boost vector for transforms between LAB <-> NRF (nucleon rest frame)
82  GHepParticle * nuc = evrec->HitNucleon();
83  const TLorentzVector & p4nuc_lab = *(nuc->P4());
84  TVector3 beta = p4nuc_lab.BoostVector();
85 
86  // Get the neutrino 4-momentum at the LAB
87  GHepParticle * nu = evrec->Probe();
88  const TLorentzVector & p4v_lab = *(nu->P4());
89 
90  // Get the neutrino 4-momentum at the NRF
91  TLorentzVector p4v_nrf = p4v_lab;
92  p4v_nrf.Boost(-1.*beta);
93 
94  // Generate the photon cos(theta) with respect to the neutrino direction
95  // (at NRF) and a uniform azimuthal angle phi
96  double costheta_gamma = -1.0 + 2.0 * rnd->RndKine().Rndm();
97  double phi_gamma = 2.0 * kPi * rnd->RndKine().Rndm();
98 
99  // Generate the fraction of the neutrino energy taken by the photon
100  double efrac_gamma = rnd->RndKine().Rndm();
101 
102  // Calculate the photon energy at the NRF
103  double Ev_nrf = p4v_nrf.Energy();
104  double Egamma_nrf = Ev_nrf * efrac_gamma;
105 
106  // Calculate the photon momentum components at a rotated NRF (NRF')
107  // where z is along the neutrino direction
108  double sintheta_gamma = TMath::Sqrt(1-TMath::Power(costheta_gamma,2));
109  double pgamma_nrf_p = Egamma_nrf * costheta_gamma; // p(//)
110  double pgamma_nrf_t = Egamma_nrf * sintheta_gamma; // p(-|)
111  double px = pgamma_nrf_t * TMath::Sin(phi_gamma);
112  double py = pgamma_nrf_t * TMath::Cos(phi_gamma);
113  double pz = pgamma_nrf_p;
114 
115  // Take a unit vector along the neutrino direction @ the NRF
116  TVector3 unit_nudir = p4v_nrf.Vect().Unit();
117 
118  // Rotate the photon momentum vector from NRF' -> NRF
119  TVector3 p3gamma_nrf(px,py,pz);
120  p3gamma_nrf.RotateUz(unit_nudir);
121 
122  // Get the photon 4-momentum back at the LAB
123  TLorentzVector p4gamma_lab(p3gamma_nrf, Egamma_nrf);
124  p4gamma_lab.Boost(beta);
125 
126  // Add the photon at the event record
127  const TLorentzVector & vtx = *(nu->X4());
128  GHepParticle p(kPdgGamma,kIStStableFinalState,0,-1,-1,-1,p4gamma_lab,vtx);
129  evrec->AddParticle(p);
130 }
131 //___________________________________________________________________________
133 {
134 // Adding the final state neutrino
135 // Just use 4-momentum conservation (init_neutrino = photon + final_neutrino)
136 
137  LOG("AMNuGammaGenerator", pINFO) << "Adding final state neutrino";
138 
139  GHepParticle * nu = evrec->Probe(); // incoming v
140  GHepParticle * gamma = evrec->Particle(nu->FirstDaughter()); // gamma
141  assert(nu);
142  assert(gamma);
143 
144  const TLorentzVector & vtx = *(nu->X4()); // vtx
145 
146  const TLorentzVector & p4nu_lab = *(nu->P4());
147  const TLorentzVector & p4gamma_lab = *(gamma->P4());
148  TLorentzVector p4_lab = p4nu_lab - p4gamma_lab;
149 
150  GHepParticle p(nu->Pdg(), kIStStableFinalState, 0,-1,-1,-1, p4_lab, vtx);
151  evrec->AddParticle(p);
152 }
153 //___________________________________________________________________________
155 {
156 // Adding the recoil nucleon.
157 // Recoil is negligible at the H3 model. However, for nuclear targets, the
158 // hit nucleon was off-the-mass-shell (bound). Doing some magic here to bring
159 // it back on-the-mass-shell so that it can be INTRANUKE'ed and appear in
160 // the final state. The nucleon will keep its original (Fermi) 3-momentum but
161 // take some energy back from the remnant nucleus. No such tweaking takes
162 // place for free nucleon targets.
163 
164  LOG("AMNuGammaGenerator", pINFO) << "Adding recoil nucleon";
165 
166  // Get the interaction target
167  GHepParticle * tgt_nucleus = evrec->TargetNucleus();
168  bool is_nuclear_target = (tgt_nucleus!=0);
169 
170  // Get the hit nucleon
171  GHepParticle * hitnuc = evrec->HitNucleon();
172  assert(hitnuc);
173 
174  // Get the hit nucleon pdg code (= recoil nucleon pdg code)
175  int pdgc = hitnuc->Pdg();
176 
177  // Get the hit nucleon 4-momentum (LAB)
178  const TLorentzVector & p4n = *(hitnuc->P4());
179  TLorentzVector p4(p4n);
180 
181  // Tweak the 4-momentum to bring the recoil nucleon on-the-mass-shell
182  // (for nuclear targets only)
183  if (is_nuclear_target) {
184  double p = p4.Vect().Mag();
185  double m = PDGLibrary::Instance()->Find(pdgc)->Mass();
186  double E = TMath::Sqrt(m*m+p*p);
187  p4.SetE(E);
188  }
189 
190  // Get the vtx position
191  GHepParticle * neutrino = evrec->Probe();
192  const TLorentzVector & vtx = *(neutrino->X4());
193 
194  // Add the recoil nucleon at the event record
195  GHepStatus_t ist = (is_nuclear_target) ?
197  int mom = evrec->HitNucleonPosition();
198 
199  LOG("AMNuGammaGenerator", pINFO)
200  << "Adding recoil baryon [pdgc = " << pdgc << "]";
201 
202  GHepParticle p(pdgc, ist, mom,-1,-1,-1, p4, vtx);
203 // double w = hitnuc->RemovalEnergy();
204 /// p.SetRemovalEnergy(w);
205  evrec->AddParticle(p);
206 }
207 //___________________________________________________________________________
209 {
210 // Add the remnant nuclear target at the GHEP record
211 
212  LOG("AMNuGammaGenerator", pINFO) << "Adding final state nucleus";
213 
214  // Skip for non nuclear targets
215  GHepParticle * nucleus = evrec->TargetNucleus();
216  if (!nucleus) {
217  LOG("AMNuGammaGenerator", pDEBUG)
218  << "No nucleus in the initial state - no remnant nucleus to add in the f/s";
219  return;
220  }
221 
222  // Compute A,Z for final state nucleus & get its PDG code and its mass
223  //GHepParticle * nucleon = evrec->HitNucleon();
224 
225  GHepParticle * hit_nucleon = evrec->HitNucleon(); // hit
226  GHepParticle * rec_nucleon = evrec->Particle(hit_nucleon->FirstDaughter()); // recoil
227 
228  assert(rec_nucleon);
229  int npdgc = rec_nucleon->Pdg();
230  bool is_p = pdg::IsProton(npdgc);
231  int A = nucleus->A();
232  int Z = nucleus->Z();
233  if (is_p) Z--;
234  A--;
235  int ipdgc = pdg::IonPdgCode(A, Z);
236  TParticlePDG * remnant = PDGLibrary::Instance()->Find(ipdgc);
237  if(!remnant) {
238  LOG("AMNuGammaGenerator", pFATAL)
239  << "No particle with [A = " << A << ", Z = " << Z
240  << ", pdgc = " << ipdgc << "] in PDGLibrary!";
241  assert(remnant);
242  }
243 
244  // Figure out the remnant 4-momentum
245  double Mf = remnant->Mass();
246  double Mf2 = TMath::Power(Mf,2);
247  double px = -1.* rec_nucleon->Px();
248  double py = -1.* rec_nucleon->Py();
249  double pz = -1.* rec_nucleon->Pz();
250  double E = TMath::Sqrt(Mf2 + rec_nucleon->P4()->Vect().Mag2());
251  E += (hit_nucleon->P4()->E() - rec_nucleon->P4()->E());
252 
253  //-- Add the nucleus to the event record
254  LOG("AMNuGammaGenerator", pINFO)
255  << "Adding nucleus [A = " << A << ", Z = " << Z
256  << ", pdgc = " << ipdgc << "]";
257  int mom = evrec->TargetNucleusPosition();
258  evrec->AddParticle(
259  ipdgc,kIStStableFinalState, mom,-1,-1,-1, px,py,pz,E, 0,0,0,0);
260 }
261 //___________________________________________________________________________
const double kPi
int Z(void) const
Basic constants.
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:148
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
int FirstDaughter(void) const
Definition: GHepParticle.h:69
enum genie::EGHepStatus GHepStatus_t
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
virtual int HitNucleonPosition(void) const
Definition: GHepRecord.cxx:454
Definition: config.py:1
Double_t beta
void AddTargetRemnant(GHepRecord *event_rec) const
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
void AddPhoton(GHepRecord *event_rec) const
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
Float_t Z
Definition: plot.C:38
int Pdg(void) const
Definition: GHepParticle.h:64
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Float_t E
Definition: plot.C:20
const int kPdgGamma
Definition: PDGCodes.h:166
void AddFinalStateNeutrino(GHepRecord *event_rec) const
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
#define pINFO
Definition: Messenger.h:63
TRandom3 & RndKine(void) const
rnd number generator used by kinematics generators
Definition: RandomGen.h:51
void AddRecoilNucleon(GHepRecord *event_rec) const
static const double A
Definition: Units.h:82
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:350
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
assert(nhit_max >=nhit_nbins)
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
int A(void) const
void ProcessEventRecord(GHepRecord *event_rec) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
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:40
virtual int TargetNucleusPosition(void) const
Definition: GHepRecord.cxx:406
#define pDEBUG
Definition: Messenger.h:64
double Py(void) const
Get Py.
Definition: GHepParticle.h:90