Generates nuclear de-excitation gamma rays. More...

#include <NucDeExcitationSim.h>

Inheritance diagram for genie::NucDeExcitationSim:

Public Member Functions
	NucDeExcitationSim ()

	NucDeExcitationSim (string config)

	~NucDeExcitationSim ()

void	ProcessEventRecord (GHepRecord *evrec) const

Public Member Functions inherited from genie::EventRecordVisitorI
virtual	~EventRecordVisitorI ()

Public Member Functions inherited from genie::Algorithm
virtual	~Algorithm ()

virtual void	Configure (const Registry &config)

virtual void	Configure (string config)

virtual void	FindConfig (void)

virtual const Registry &	GetConfig (void) const

Registry *	GetOwnedConfig (void)

virtual const AlgId &	Id (void) const
	Get algorithm ID. More...

virtual AlgStatus_t	GetStatus (void) const
	Get algorithm status. More...

virtual bool	AllowReconfig (void) const

virtual AlgCmp_t	Compare (const Algorithm *alg) const
	Compare with input algorithm. More...

virtual void	SetId (const AlgId &id)
	Set algorithm ID. More...

virtual void	SetId (string name, string config)

const Algorithm *	SubAlg (const RgKey &registry_key) const

void	AdoptConfig (void)

void	AdoptSubstructure (void)

virtual void	Print (ostream &stream) const
	Print algorithm info. More...

Private Member Functions
void	OxygenTargetSim (GHepRecord *evrec) const

void	AddPhoton (GHepRecord *evrec, double E0, double t) const

double	PhotonEnergySmearing (double E0, double t) const

TLorentzVector	Photon4P (double E) const

Additional Inherited Members
Static Public Member Functions inherited from genie::Algorithm
static string	BuildParamVectKey (const std::string &comm_name, unsigned int i)

static string	BuildParamVectSizeKey (const std::string &comm_name)

Protected Member Functions inherited from genie::EventRecordVisitorI
	EventRecordVisitorI ()

	EventRecordVisitorI (string name)

	EventRecordVisitorI (string name, string config)

Protected Member Functions inherited from genie::Algorithm
	Algorithm ()

	Algorithm (string name)

	Algorithm (string name, string config)

void	Initialize (void)

void	DeleteConfig (void)

void	DeleteSubstructure (void)

Registry *	ExtractLocalConfig (const Registry &in) const

Registry *	ExtractLowerConfig (const Registry &in, const string &alg_key) const
	Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...

template<class T >
bool	GetParam (const RgKey &name, T &p, bool is_top_call=true) const

template<class T >
bool	GetParamDef (const RgKey &name, T &p, const T &def) const

template<class T >
int	GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
	Handle to load vectors of parameters. More...

int	GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const

int	AddTopRegistry (Registry *rp, bool owns=true)
	add registry with top priority, also update ownership More...

int	AddLowRegistry (Registry *rp, bool owns=true)
	add registry with lowest priority, also update ownership More...

int	MergeTopRegistry (const Registry &r)

int	AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
	Add registries with top priority, also udated Ownerships. More...

Protected Attributes inherited from genie::Algorithm
bool	fAllowReconfig

bool	fOwnsSubstruc
	true if it owns its substructure (sub-algs,...) More...

AlgId	fID
	algorithm name and configuration set More...

vector< Registry * >	fConfVect

vector< bool >	fOwnerships
	ownership for every registry in fConfVect More...

AlgStatus_t	fStatus
	algorithm execution status More...

AlgMap *	fOwnedSubAlgMp
	local pool for owned sub-algs (taken out of the factory pool) More...

Detailed Description

Generates nuclear de-excitation gamma rays.

Author: Costas Andreopoulos <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

16O: H.Ejiri,Phys.Rev.C48,1442(1993); K.Kobayashi et al., Nucl.Phys.B (proc Suppl) 139 (2005)

March 05, 2008

Definition at line 31 of file NucDeExcitationSim.h.

Constructor & Destructor Documentation

NucDeExcitationSim::NucDeExcitationSim ( )

Definition at line 43 of file NucDeExcitationSim.cxx.

                                        :
 EventRecordVisitorI("genie::NucDeExcitationSim")
 {
 
 }

NucDeExcitationSim::NucDeExcitationSim ( string config )

Definition at line 49 of file NucDeExcitationSim.cxx.

                                                     :
 EventRecordVisitorI("genie::NucDeExcitationSim", config)
 {
 
 }

NucDeExcitationSim::~NucDeExcitationSim ( )

Definition at line 55 of file NucDeExcitationSim.cxx.

 {
 
 }

Member Function Documentation

void NucDeExcitationSim::AddPhoton	(	GHepRecord *	evrec,
		double	E0,
		double	t
	)		const

private

Definition at line 309 of file NucDeExcitationSim.cxx.

 {
 // Add a photon at the event record & recoil the remnant nucleus so as to
 // conserve energy/momenta
 //
   double E = (dt>0) ? this->PhotonEnergySmearing(E0, dt) : E0;
 
   LOG("NucDeEx", pNOTICE)
     << "Adding a " << E/units::MeV << " MeV photon from nucl. deexcitation";
 
   GHepParticle * target  = evrec->Particle(1);
   GHepParticle * remnant = 0;
   for(int i = target->FirstDaughter(); i <= target->LastDaughter(); i++) {
     remnant  = evrec->Particle(i);
     if(pdg::IsIon(remnant->Pdg())) break;
   }
 
   TLorentzVector x4(0,0,0,0);
   TLorentzVector p4 = this->Photon4P(E);
   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??)
   evrec->AddParticle(gamma);
 
 
   remnant->SetPx     ( remnant->Px() - p4.Px() );
   remnant->SetPy     ( remnant->Py() - p4.Py() );
   remnant->SetPz     ( remnant->Pz() - p4.Pz() );
   remnant->SetEnergy ( remnant->E()  - p4.E()  );
 }

void NucDeExcitationSim::OxygenTargetSim ( GHepRecord * evrec ) const

private

Definition at line 78 of file NucDeExcitationSim.cxx.

 {
   LOG("NucDeEx", pNOTICE)
      << "Simulating nuclear de-excitation gamma rays for Oxygen target";
 
   //LOG("NucDeEx", pNOTICE) << *evrec;
 
   GHepParticle * hitnuc = evrec->HitNucleon();
   if(!hitnuc) return;
 
   bool p_hole = (hitnuc->Pdg() == kPdgProton);
   double dt   = -1;
 
   RandomGen * rnd = RandomGen::Instance();
 
   //
   // ****** P-Hole
   //
   if (p_hole) {
     //
     // * Define all the data required for simulating deexcitations of p-hole states
     //
 
     // > probabilities for creating a p-hole in the P1/2, P3/2, S1/2 shells
     double Pp12 = 0.25;              // P1/2
     double Pp32 = 0.47;              // P3/2
     double Ps12 = 1. - Pp12 - Pp32;  // S1/2
 
     // > excited state energy levels & probabilities for P3/2-shell p-holes
     const int np32 = 3;
     double p32Elv[np32] = { 0.00632, 0.00993, 0.01070 };
     double p32Plv[np32] = { 0.872,   0.064,   0.064   };
     // - probabilities for deexcitation modes of P3/2-shell p-hole state '1'
     double p32Plv1_1gamma  = 0.78;  // prob to decay via 1 gamma
     double p32Plv1_cascade = 0.22;  // prob to decay via gamma cascade
 
     // > excited state energy levels & probabilities for S1/2-shell p-holes
     const int ns12 = 11;
     double s12Elv[ns12] = {
                0.00309, 0.00368, 0.00385, 0.00444, 0.00492,
                0.00511, 0.00609, 0.00673, 0.00701, 0.00703, 0.00734 };
     double s12Plv[ns12] = {
                0.0625,  0.1875,  0.075,   0.1375,  0.1375,
                0.0125,  0.0125,  0.075,   0.0563,  0.0563,  0.1874  };
     // - gamma energies and probabilities for S1/2-shell p-hole excited
     //   states '2','7' and '10' with >1 deexcitation modes
     const int ns12lv2 = 3;
     double s12Elv2[ns12lv2]    = { 0.00309, 0.00369, 0.00385 };
     double s12Plv2[ns12lv2]    = { 0.013,   0.360,   0.625   };
     const int ns12lv7 = 2;
     double s12Elv7[ns12lv7]    = { 0.00609, 0.00673 };
     double s12Plv7[ns12lv7]    = { 0.04,    0.96    };
     const int ns12lv10 = 3;
     double s12Elv10[ns12lv10]  = { 0.00609, 0.00673, 0.00734 };
     double s12Plv10[ns12lv10]  = { 0.050,   0.033,   0.017   };
 
     // Select one of the P1/2, P3/2 or S1/2
     double rshell = rnd->RndDec().Rndm();
     //
     // >> P1/2 shell
     //
     if(rshell < Pp12) {
         LOG("NucDeEx", pNOTICE)
           << "Hit nucleon left a P1/2 shell p-hole. Remnant is at g.s.";
         return;
     }
     //
     // >> P3/2 shell
     //
     else
     if(rshell < Pp12 + Pp32) {
         LOG("NucDeEx", pNOTICE)
             << "Hit nucleon left a P3/2 shell p-hole";
         // Select one of the excited states
         double rdecmode  = rnd->RndDec().Rndm();
         double prob_sum  = 0;
         int    sel_state = -1;
         for(int istate=0; istate<np32; istate++) {
             prob_sum += p32Plv[istate];
             if(rdecmode < prob_sum) {
               sel_state = istate;
               break;
             }
         }
         LOG("NucDeEx", pNOTICE)
             << "Selected P3/2 excited state = " << sel_state;
 
         // Decay that excited state
         // >> 6.32 MeV state
         if(sel_state==0) {
             this->AddPhoton(evrec, p32Elv[0], dt);
         }
         // >> 9.93 MeV state
         else
         if(sel_state==1) {
             double r = rnd->RndDec().Rndm();
             // >>> emit a single gamma
             if(r < p32Plv1_1gamma) {
                this->AddPhoton(evrec, p32Elv[1], dt);
             }
             // >>> emit a cascade of gammas
             else
             if(r < p32Plv1_1gamma + p32Plv1_cascade) {
                this->AddPhoton(evrec, p32Elv[1],           dt);
                this->AddPhoton(evrec, p32Elv[1]-p32Elv[0], dt);
             }
         }
         // >> 10.7 MeV state
         else
         if(sel_state==2) {
            // Above the particle production threshold - need to emit
            // a 0.5 MeV kinetic energy proton.
            // Will neglect that given that it is a very low energy
            // kinetic energy nucleon and the intranuke break-up nucleon
            // cross sections are already tuned.
            return;
         }
     } //p3/2
     //
     // >> S1/2 shell
     //
     else if (rshell < Pp12 + Pp32 + Ps12) {
         LOG("NucDeEx", pNOTICE)
             << "Hit nucleon left an S1/2 shell p-hole";
         // Select one of the excited states caused by a S1/2 shell hole
         double rdecmode  = rnd->RndDec().Rndm();
         double prob_sum  = 0;
         int    sel_state = -1;
         for(int istate=0; istate<ns12; istate++) {
             prob_sum += s12Plv[istate];
             if(rdecmode < prob_sum) {
               sel_state = istate;
               break;
             }
         }
         LOG("NucDeEx", pNOTICE)
             << "Selected S1/2 excited state = " << sel_state;
 
         // Decay that excited state
         bool multiple_decay_modes =
               (sel_state==2 || sel_state==7 || sel_state==10);
         if(!multiple_decay_modes) {
           this->AddPhoton(evrec, s12Elv[sel_state], dt);
         } else {
           int ndec = -1;
           double * pdec = 0, * edec = 0;
           switch(sel_state) {
            case(2) :
               ndec = ns12lv2;  pdec = s12Plv2;  edec = s12Elv2;
               break;
            case(7) :
               ndec = ns12lv7;  pdec = s12Plv7;  edec = s12Elv7;
               break;
            case(10) :
               ndec = ns12lv10; pdec = s12Plv10; edec = s12Elv10;
               break;
            default:
              return;
           }
           double r = rnd->RndDec().Rndm();
           double decmode_prob_sum = 0;
           int sel_decmode = -1;
           for(int idecmode=0; idecmode < ndec; idecmode++) {
              decmode_prob_sum += pdec[idecmode];
              if(r < decmode_prob_sum) {
                  sel_decmode = idecmode;
                  break;
              }
           }
           if(sel_decmode == -1) return;
           this->AddPhoton(evrec, edec[sel_decmode], dt);
         }//mult.dec.ch
 
     } // s1/2
     else {
     }
   } // p-hole
 
   //
   // ****** n-hole
   //
   else {
     //
     // * Define all the data required for simulating deexcitations of n-hole states
     //
 
     // > probabilities for creating a n-hole in the P1/2, P3/2, S1/2 shells
     double Pp12 = 0.25;  // P1/2
     double Pp32 = 0.44;  // P3/2
     double Ps12 = 0.09;  // S1/2
     //>
     double p32Elv = 0.00618;
     //>
     double s12Elv = 0.00703;
     double s12Plv = 0.222;
 
     // Select one of the P1/2, P3/2 or S1/2
     double rshell = rnd->RndDec().Rndm();
     //
     // >> P1/2 shell
     //
     if(rshell < Pp12) {
         LOG("NucDeEx", pNOTICE)
           << "Hit nucleon left a P1/2 shell n-hole. Remnant is at g.s.";
         return;
     }
     //
     // >> P3/2 shell
     //
     else
     if(rshell < Pp12 + Pp32) {
         LOG("NucDeEx", pNOTICE)
             << "Hit nucleon left a P3/2 shell n-hole";
         this->AddPhoton(evrec, p32Elv, dt);
     }
     //
     // >> S1/2 shell
     //
     else
     if(rshell < Pp12 + Pp32 + Ps12) {
         LOG("NucDeEx", pNOTICE)
             << "Hit nucleon left a S1/2 shell n-hole";
         // only one of the deexcitation modes involve a (7.03 MeV) photon
         double r = rnd->RndDec().Rndm();
         if(r < s12Plv) this->AddPhoton(evrec, s12Elv,dt);
     }
     else {
     }
   } //n-hole
 }

TLorentzVector NucDeExcitationSim::Photon4P ( double E ) const

private

Definition at line 356 of file NucDeExcitationSim.cxx.

 {
 // Generate a photon 4p
 
   RandomGen * rnd = RandomGen::Instance();
 
   double costheta = -1. + 2. * rnd->RndDec().Rndm();
   double sintheta = TMath::Sqrt(TMath::Max(0., 1.-TMath::Power(costheta,2)));
   double phi      = 2*kPi * rnd->RndDec().Rndm();
   double cosphi   = TMath::Cos(phi);
   double sinphi   = TMath::Sin(phi);
 
   double px = E * sintheta * cosphi;
   double py = E * sintheta * sinphi;
   double pz = E * costheta;
 
   TLorentzVector p4(px,py,pz,E);
   return p4;
 }

double NucDeExcitationSim::PhotonEnergySmearing	(	double	E0,
		double	t
	)		const

private

Definition at line 339 of file NucDeExcitationSim.cxx.

 {
 // Returns the smeared energy of the emitted gamma
 // E0 : energy of the excited state (GeV)
 // dt: excited state lifetime (sec)
 //
   double dE = kPlankConstant / (dt*units::s);
 
   RandomGen * rnd = RandomGen::Instance();
   double E = rnd->RndDec().Gaus(E0 /*mean*/, dE /*sigma*/);
 
   LOG("NucDeEx", pNOTICE)
      << "<E> = " << E0 << ", dE = " << dE << " -> E = " << E;
 
   return E;
 }

void NucDeExcitationSim::ProcessEventRecord ( GHepRecord * evrec ) const

virtual

Implements genie::EventRecordVisitorI.

Definition at line 60 of file NucDeExcitationSim.cxx.

 {
   LOG("NucDeEx", pNOTICE)
      << "Simulating nuclear de-excitation gamma rays";
 
   GHepParticle * nucltgt = evrec->TargetNucleus();
   if (!nucltgt) {
     LOG("NucDeEx", pINFO)
       << "No nuclear target found - Won't simulate nuclear de-excitation";
     return;
   }
 
   if(nucltgt->Z()==8) this->OxygenTargetSim(evrec);
 
   LOG("NucDeEx", pINFO)
      << "Done with this event";
 }

The documentation for this class was generated from the following files:

Generator/src/Physics/NuclearDeExcitation/NucDeExcitationSim.h
Generator/src/Physics/NuclearDeExcitation/NucDeExcitationSim.cxx

Public Member Functions

Private Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation