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