NucDeExcitationSim.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2020, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
7  University of Liverpool & STFC Rutherford Appleton Laboratory
8 */
9 //____________________________________________________________________________
10 
11 #include <cstdlib>
12 #include <sstream>
13 
14 #include <TMath.h>
15 
17 #include "Framework/Conventions/GBuild.h"
34 
35 using std::ostringstream;
36 
37 using namespace genie;
38 using namespace genie::utils;
39 using namespace genie::constants;
40 using namespace genie::controls;
41 
42 //___________________________________________________________________________
44 EventRecordVisitorI("genie::NucDeExcitationSim")
45 {
46 
47 }
48 //___________________________________________________________________________
50 EventRecordVisitorI("genie::NucDeExcitationSim", config)
51 {
52 
53 }
54 //___________________________________________________________________________
56 {
57 
58 }
59 //___________________________________________________________________________
61 {
62  LOG("NucDeEx", pNOTICE)
63  << "Simulating nuclear de-excitation gamma rays";
64 
65  GHepParticle * nucltgt = evrec->TargetNucleus();
66  if (!nucltgt) {
67  LOG("NucDeEx", pINFO)
68  << "No nuclear target found - Won't simulate nuclear de-excitation";
69  return;
70  }
71 
72  if(nucltgt->Z()==8) this->OxygenTargetSim(evrec);
73 
74  LOG("NucDeEx", pINFO)
75  << "Done with this event";
76 }
77 //___________________________________________________________________________
79 {
80  LOG("NucDeEx", pNOTICE)
81  << "Simulating nuclear de-excitation gamma rays for Oxygen target";
82 
83  //LOG("NucDeEx", pNOTICE) << *evrec;
84 
85  GHepParticle * hitnuc = evrec->HitNucleon();
86  if(!hitnuc) return;
87 
88  bool p_hole = (hitnuc->Pdg() == kPdgProton);
89  double dt = -1;
90 
92 
93  //
94  // ****** P-Hole
95  //
96  if (p_hole) {
97  //
98  // * Define all the data required for simulating deexcitations of p-hole states
99  //
100 
101  // > probabilities for creating a p-hole in the P1/2, P3/2, S1/2 shells
102  double Pp12 = 0.25; // P1/2
103  double Pp32 = 0.47; // P3/2
104  double Ps12 = 1. - Pp12 - Pp32; // S1/2
105 
106  // > excited state energy levels & probabilities for P3/2-shell p-holes
107  const int np32 = 3;
108  double p32Elv[np32] = { 0.00632, 0.00993, 0.01070 };
109  double p32Plv[np32] = { 0.872, 0.064, 0.064 };
110  // - probabilities for deexcitation modes of P3/2-shell p-hole state '1'
111  double p32Plv1_1gamma = 0.78; // prob to decay via 1 gamma
112  double p32Plv1_cascade = 0.22; // prob to decay via gamma cascade
113 
114  // > excited state energy levels & probabilities for S1/2-shell p-holes
115  const int ns12 = 11;
116  double s12Elv[ns12] = {
117  0.00309, 0.00368, 0.00385, 0.00444, 0.00492,
118  0.00511, 0.00609, 0.00673, 0.00701, 0.00703, 0.00734 };
119  double s12Plv[ns12] = {
120  0.0625, 0.1875, 0.075, 0.1375, 0.1375,
121  0.0125, 0.0125, 0.075, 0.0563, 0.0563, 0.1874 };
122  // - gamma energies and probabilities for S1/2-shell p-hole excited
123  // states '2','7' and '10' with >1 deexcitation modes
124  const int ns12lv2 = 3;
125  double s12Elv2[ns12lv2] = { 0.00309, 0.00369, 0.00385 };
126  double s12Plv2[ns12lv2] = { 0.013, 0.360, 0.625 };
127  const int ns12lv7 = 2;
128  double s12Elv7[ns12lv7] = { 0.00609, 0.00673 };
129  double s12Plv7[ns12lv7] = { 0.04, 0.96 };
130  const int ns12lv10 = 3;
131  double s12Elv10[ns12lv10] = { 0.00609, 0.00673, 0.00734 };
132  double s12Plv10[ns12lv10] = { 0.050, 0.033, 0.017 };
133 
134  // Select one of the P1/2, P3/2 or S1/2
135  double rshell = rnd->RndDec().Rndm();
136  //
137  // >> P1/2 shell
138  //
139  if(rshell < Pp12) {
140  LOG("NucDeEx", pNOTICE)
141  << "Hit nucleon left a P1/2 shell p-hole. Remnant is at g.s.";
142  return;
143  }
144  //
145  // >> P3/2 shell
146  //
147  else
148  if(rshell < Pp12 + Pp32) {
149  LOG("NucDeEx", pNOTICE)
150  << "Hit nucleon left a P3/2 shell p-hole";
151  // Select one of the excited states
152  double rdecmode = rnd->RndDec().Rndm();
153  double prob_sum = 0;
154  int sel_state = -1;
155  for(int istate=0; istate<np32; istate++) {
156  prob_sum += p32Plv[istate];
157  if(rdecmode < prob_sum) {
158  sel_state = istate;
159  break;
160  }
161  }
162  LOG("NucDeEx", pNOTICE)
163  << "Selected P3/2 excited state = " << sel_state;
164 
165  // Decay that excited state
166  // >> 6.32 MeV state
167  if(sel_state==0) {
168  this->AddPhoton(evrec, p32Elv[0], dt);
169  }
170  // >> 9.93 MeV state
171  else
172  if(sel_state==1) {
173  double r = rnd->RndDec().Rndm();
174  // >>> emit a single gamma
175  if(r < p32Plv1_1gamma) {
176  this->AddPhoton(evrec, p32Elv[1], dt);
177  }
178  // >>> emit a cascade of gammas
179  else
180  if(r < p32Plv1_1gamma + p32Plv1_cascade) {
181  this->AddPhoton(evrec, p32Elv[1], dt);
182  this->AddPhoton(evrec, p32Elv[1]-p32Elv[0], dt);
183  }
184  }
185  // >> 10.7 MeV state
186  else
187  if(sel_state==2) {
188  // Above the particle production threshold - need to emit
189  // a 0.5 MeV kinetic energy proton.
190  // Will neglect that given that it is a very low energy
191  // kinetic energy nucleon and the intranuke break-up nucleon
192  // cross sections are already tuned.
193  return;
194  }
195  } //p3/2
196  //
197  // >> S1/2 shell
198  //
199  else if (rshell < Pp12 + Pp32 + Ps12) {
200  LOG("NucDeEx", pNOTICE)
201  << "Hit nucleon left an S1/2 shell p-hole";
202  // Select one of the excited states caused by a S1/2 shell hole
203  double rdecmode = rnd->RndDec().Rndm();
204  double prob_sum = 0;
205  int sel_state = -1;
206  for(int istate=0; istate<ns12; istate++) {
207  prob_sum += s12Plv[istate];
208  if(rdecmode < prob_sum) {
209  sel_state = istate;
210  break;
211  }
212  }
213  LOG("NucDeEx", pNOTICE)
214  << "Selected S1/2 excited state = " << sel_state;
215 
216  // Decay that excited state
217  bool multiple_decay_modes =
218  (sel_state==2 || sel_state==7 || sel_state==10);
219  if(!multiple_decay_modes) {
220  this->AddPhoton(evrec, s12Elv[sel_state], dt);
221  } else {
222  int ndec = -1;
223  double * pdec = 0, * edec = 0;
224  switch(sel_state) {
225  case(2) :
226  ndec = ns12lv2; pdec = s12Plv2; edec = s12Elv2;
227  break;
228  case(7) :
229  ndec = ns12lv7; pdec = s12Plv7; edec = s12Elv7;
230  break;
231  case(10) :
232  ndec = ns12lv10; pdec = s12Plv10; edec = s12Elv10;
233  break;
234  default:
235  return;
236  }
237  double r = rnd->RndDec().Rndm();
238  double decmode_prob_sum = 0;
239  int sel_decmode = -1;
240  for(int idecmode=0; idecmode < ndec; idecmode++) {
241  decmode_prob_sum += pdec[idecmode];
242  if(r < decmode_prob_sum) {
243  sel_decmode = idecmode;
244  break;
245  }
246  }
247  if(sel_decmode == -1) return;
248  this->AddPhoton(evrec, edec[sel_decmode], dt);
249  }//mult.dec.ch
250 
251  } // s1/2
252  else {
253  }
254  } // p-hole
255 
256  //
257  // ****** n-hole
258  //
259  else {
260  //
261  // * Define all the data required for simulating deexcitations of n-hole states
262  //
263 
264  // > probabilities for creating a n-hole in the P1/2, P3/2, S1/2 shells
265  double Pp12 = 0.25; // P1/2
266  double Pp32 = 0.44; // P3/2
267  double Ps12 = 0.09; // S1/2
268  //>
269  double p32Elv = 0.00618;
270  //>
271  double s12Elv = 0.00703;
272  double s12Plv = 0.222;
273 
274  // Select one of the P1/2, P3/2 or S1/2
275  double rshell = rnd->RndDec().Rndm();
276  //
277  // >> P1/2 shell
278  //
279  if(rshell < Pp12) {
280  LOG("NucDeEx", pNOTICE)
281  << "Hit nucleon left a P1/2 shell n-hole. Remnant is at g.s.";
282  return;
283  }
284  //
285  // >> P3/2 shell
286  //
287  else
288  if(rshell < Pp12 + Pp32) {
289  LOG("NucDeEx", pNOTICE)
290  << "Hit nucleon left a P3/2 shell n-hole";
291  this->AddPhoton(evrec, p32Elv, dt);
292  }
293  //
294  // >> S1/2 shell
295  //
296  else
297  if(rshell < Pp12 + Pp32 + Ps12) {
298  LOG("NucDeEx", pNOTICE)
299  << "Hit nucleon left a S1/2 shell n-hole";
300  // only one of the deexcitation modes involve a (7.03 MeV) photon
301  double r = rnd->RndDec().Rndm();
302  if(r < s12Plv) this->AddPhoton(evrec, s12Elv,dt);
303  }
304  else {
305  }
306  } //n-hole
307 }
308 //___________________________________________________________________________
310  GHepRecord * evrec, double E0, double dt) const
311 {
312 // Add a photon at the event record & recoil the remnant nucleus so as to
313 // conserve energy/momenta
314 //
315  double E = (dt>0) ? this->PhotonEnergySmearing(E0, dt) : E0;
316 
317  LOG("NucDeEx", pNOTICE)
318  << "Adding a " << E/units::MeV << " MeV photon from nucl. deexcitation";
319 
320  GHepParticle * target = evrec->Particle(1);
321  GHepParticle * remnant = 0;
322  for(int i = target->FirstDaughter(); i <= target->LastDaughter(); i++) {
323  remnant = evrec->Particle(i);
324  if(pdg::IsIon(remnant->Pdg())) break;
325  }
326 
327  TLorentzVector x4(0,0,0,0);
328  TLorentzVector p4 = this->Photon4P(E);
329  GHepParticle gamma(kPdgGamma, kIStStableFinalState,1,-1,-1,-1, p4, x4); // note that this assigns the parent of the photon as the initial-state nucleon/nucleus. (do we want that??)
330  evrec->AddParticle(gamma);
331 
332 
333  remnant->SetPx ( remnant->Px() - p4.Px() );
334  remnant->SetPy ( remnant->Py() - p4.Py() );
335  remnant->SetPz ( remnant->Pz() - p4.Pz() );
336  remnant->SetEnergy ( remnant->E() - p4.E() );
337 }
338 //___________________________________________________________________________
339 double NucDeExcitationSim::PhotonEnergySmearing(double E0, double dt) const
340 {
341 // Returns the smeared energy of the emitted gamma
342 // E0 : energy of the excited state (GeV)
343 // dt: excited state lifetime (sec)
344 //
345  double dE = kPlankConstant / (dt*units::s);
346 
347  RandomGen * rnd = RandomGen::Instance();
348  double E = rnd->RndDec().Gaus(E0 /*mean*/, dE /*sigma*/);
349 
350  LOG("NucDeEx", pNOTICE)
351  << "<E> = " << E0 << ", dE = " << dE << " -> E = " << E;
352 
353  return E;
354 }
355 //___________________________________________________________________________
356 TLorentzVector NucDeExcitationSim::Photon4P(double E) const
357 {
358 // Generate a photon 4p
359 
360  RandomGen * rnd = RandomGen::Instance();
361 
362  double costheta = -1. + 2. * rnd->RndDec().Rndm();
363  double sintheta = TMath::Sqrt(TMath::Max(0., 1.-TMath::Power(costheta,2)));
364  double phi = 2*kPi * rnd->RndDec().Rndm();
365  double cosphi = TMath::Cos(phi);
366  double sinphi = TMath::Sin(phi);
367 
368  double px = E * sintheta * cosphi;
369  double py = E * sintheta * sinphi;
370  double pz = E * costheta;
371 
372  TLorentzVector p4(px,py,pz,E);
373  return p4;
374 }
375 //___________________________________________________________________________
int Z(void) const
Basic constants.
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:104
double PhotonEnergySmearing(double E0, double t) const
THE MAIN GENIE PROJECT NAMESPACE
Definition: AlgCmp.h:25
double E(void) const
Get energy.
Definition: GHepParticle.h:91
void SetPz(double pz)
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
int FirstDaughter(void) const
Definition: GHepParticle.h:68
static constexpr double s
Definition: Units.h:95
TLorentzVector Photon4P(double E) const
static constexpr double MeV
Definition: Units.h:129
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
void SetPx(double px)
void ProcessEventRecord(GHepRecord *evrec) const
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
int Pdg(void) const
Definition: GHepParticle.h:63
int LastDaughter(void) const
Definition: GHepParticle.h:69
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static Config * config
Definition: config.cpp:1054
const int kPdgGamma
Definition: PDGCodes.h:189
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:286
#define pINFO
Definition: Messenger.h:62
double gamma(double KE, const simb::MCParticle *part)
Misc GENIE control constants.
void AddPhoton(GHepRecord *evrec, double E0, double t) const
bool IsIon(int pdgc)
Definition: PDGUtils.cxx:39
void OxygenTargetSim(GHepRecord *evrec) const
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:306
E
Definition: 018_def.c:13
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:491
const int kPdgProton
Definition: PDGCodes.h:81
#define pNOTICE
Definition: Messenger.h:61
void SetEnergy(double E)
static const double kPlankConstant
Definition: Constants.h:32
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
static const double kPi
Definition: Constants.h:37
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
TRandom3 & RndDec(void) const
rnd number generator used by decay models
Definition: RandomGen.h:56
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
void SetPy(double py)
Root of GENIE utility namespaces.
double Py(void) const
Get Py.
Definition: GHepParticle.h:89