TParticle Class Reference

This class is how the properties of a nucleus, its spatial distribution and energy spectrum are described in DRAGON. More...

#include <particle.h>

List of all members.

Public Member Functions

 TParticle ()
 TParticle (TParticle &part)
 TParticle (int A_, int Z_, Galaxy *gal, Input *in, int uid_prec, TXSecBase *xsecmodel, TNucleiList *l)
 ~TParticle ()
TParticle & operator+= (const TParticle &addendum)
void Evolve (vector< TParticle * >, vector< TCREvolutorBasis * >, TSpallationNetwork *spnet)
void Print (fitsfile *, double)
void PrintSpectrum (fitsfile *, double)
double FindNormalization ()
 Find normalization factor, by comparing with observed proton flux.
TSpectrum * GetSpectrum ()
double GetDensity (int ir, int iz, int ip)
vector< double > & GetDensity ()
double GetDensity (int i)
int GetUid ()
const double GetLifetime () const
const int GetDaughter () const
int GetA ()
int GetZ ()
int GetIsSec ()

Protected Member Functions

vector< double > ComputeSecondarySource (vector< TParticle * >, TSpallationNetwork *spnet)
int index (int ir, int iz, int ip)
 Convert matrix to linear representation.

Protected Attributes

int A
int Z
int uid
int dimr
int dimz
int dimE
double lifetime
int daughter
int issec
Galaxy * _fGalaxy
TInelasticCrossSection * _fInXSec
TSpectrum * sp
vector< TEnergyLoss * > eloss
vector< double > density

Detailed Description

This class is how the properties of a nucleus, its spatial distribution and energy spectrum are described in DRAGON.

Author:
Luca Maccione luca.maccione@desy.de A nucleus is defined to have mass and charge, a unique ID computed from them, possibly some decay channel and daughter nuclei. Moreover, it has a distribution in physical space and in momentum which has to be found after propagation, and its propagation occurs in a given galaxy, where the nucleus suffers energy losses and inelastic scattering. Of course, nuclei are injected by sources in the ISM with some spectrum.

Definition at line 46 of file particle.h.

Constructor & Destructor Documentation

TParticle::TParticle (  )  [inline]

Default constructor.

Definition at line 49 of file particle.h.

TParticle::TParticle ( TParticle &  part  )  [inline]

Copy constructor.

Definition at line 50 of file particle.h.

References _fGalaxy, _fInXSec, density, eloss, issec, and sp.

TParticle::TParticle ( int  A_,
int  Z_,
Galaxy *  gal,
Input *  in,
int  uid_prec,
TXSecBase *  xsecmodel,
TNucleiList *  l 
)

The constructor normally used in DRAGON

Parameters:
A_ Mass number
Z_ Charge
gal A model for the galaxy
in User input

Definition at line 22 of file particle.cc.

References _fGalaxy, _fInXSec, A, daughter, density, dimE, dimr, dimz, Input::dmmode, eloss, Galaxy::GetBField(), Galaxy::GetCoordinates(), Galaxy::GetGas(), Galaxy::GetInjLowSpectrum(), Galaxy::GetInjSpectrum(), Galaxy::GetISRF(), TNucleiList::GetLifeTime(), Galaxy::GetTotalGas(), issec, lifetime, Input::mx, Myr, Input::sigmav, sp, Input::taudec, uid, and Z.

TParticle::~TParticle (  ) 

Destructor

Definition at line 62 of file particle.cc.

References _fInXSec, density, eloss, and sp.

00062                       { 
00063   if (sp) delete sp; 
00064   if (_fInXSec) delete _fInXSec;
00065   for (vector<TEnergyLoss*>::iterator i = eloss.begin(); i != eloss.end(); ++i) {
00066     if (*i) delete *i;
00067   }
00068   eloss.clear();
00069   density.clear();
00070 }

Member Function Documentation

vector< double > TParticle::ComputeSecondarySource ( vector< TParticle * >  part,
TSpallationNetwork *  spnet 
) [protected]

< Here also inverse beta-decay (K-electron capture) must be taken into account (which is not, up to now).

Parameters:
part Array of already propagated heavier nuclei, that can possibly generate present nucleus during spallation

Definition at line 115 of file particle.cc.

References _fGalaxy, A, density, dimE, dimr, dimz, TGrid::GetBeta(), Galaxy::GetCoordinates(), TGrid::GetGamma(), TGas::GetGas(), Galaxy::GetTotalGas(), TSpallationNetwork::GetXSec(), TSpallationNetwork::GetXSecApEl(), TGrid::index(), issec, and uid.

Referenced by Evolve().

00115                                                                                                    {
00118   vector<double> result(density.size(), 0.0);
00119   // If it is the only particle to be propagated, or if it is primary electrons, do not add secondary contribution
00120   if (part.size() == 1 || (uid == -1000 && !issec)) return result; 
00121 
00122   TGas* totalgas = _fGalaxy->GetTotalGas();
00123   TGrid* coord = _fGalaxy->GetCoordinates();
00124 
00125   if (issec && uid != -1000) { // Secondary protons or tertiary antiprotons: they are next to their primary brothers...
00126 
00127     vector<TParticle*>::iterator iprim = part.end()-2; // Primary species
00128 
00129     vector<double> beta(coord->GetBeta());
00130     vector<double> spall_spectrum( spnet->GetXSec(uid, uid) );
00131 
00132     for (int k = 0; k < dimr; ++k) {
00133       for (int l = 0; l < dimz; ++l) {
00134         int indspat = coord->index(k,l);
00135         double gasdensity = totalgas->GetGas(indspat);
00136 
00137         for (int i = 0; i < dimE; ++i)  {
00138           double betaigasdensity = beta[i]*gasdensity;
00139           int ind = indspat*dimE+i;       
00140           for (int ii = i+1; ii < dimE; ++ii) result[ind] += (*iprim)->GetDensity(indspat*dimE+ii) * betaigasdensity * spall_spectrum[ii];
00141         }
00142       }
00143     }
00144   } // if (issec && uid != -1000)
00145   else if (uid > 1000) { // nuclei spallation
00146 
00147     vector<double> gamma(coord->GetGamma());
00148 
00149     for (vector<TParticle*>::iterator ipart = part.begin(); ipart != part.end()-1; ++ipart) {
00150       vector<double> spall_spectrum( spnet->GetXSec( (*ipart)->GetUid(), uid ) );
00151       if (spall_spectrum.size() == dimE || (*ipart)->GetDaughter() == uid) {
00152         
00153         double Afactor = double((*ipart)->GetA())/double(A); // To convert between kinetic energy per nucleon (which is approx. constant) to momentum
00154 
00155         for (int k = 0; k < dimr; ++k) {
00156           for (int l = 0; l < dimz; ++l) {
00157             int indspat = coord->index(k,l);
00158             double Afactorgasdensity = Afactor*totalgas->GetGas(indspat);
00159             for (int i = 0; i < dimE; ++i) {
00160               int ind = indspat*dimE+i;
00161               if (spall_spectrum.size() == dimE) result[ind] += Afactorgasdensity * spall_spectrum[i] * (*ipart)->GetDensity(ind); // spallation
00162               if ((*ipart)->GetDaughter() == uid) result[ind] += (*ipart)->GetDensity(ind)/(*ipart)->GetLifetime()/gamma[i]; // decay
00163             }
00164           }
00165         }
00166       }
00167     }
00168   } // else if (uid > 1000)
00169   else { // Secondary antiprotons, electrons and positrons
00170 
00171     for (vector<TParticle*>::iterator ipart = part.begin(); ipart != part.end()-1; ++ipart) {
00172       if ((*ipart)->GetUid() > 2004 || (*ipart)->GetUid() < 0) continue;
00173       for (int i = 0; i < dimE; i++) {
00174         vector<double> spall_spectrum( spnet->GetXSecApEl((*ipart)->GetUid(), uid, i) );
00175         if (spall_spectrum.size() == dimE) {
00176 
00177           for (int k = 0; k < dimr; k++) {
00178             for (int l = 0; l < dimz; l++) {
00179               int indspat = coord->index(k,l);
00180               double gasdensity = totalgas->GetGas(indspat);
00181               int ind = indspat*dimE+i;   
00182               // ENERGY INTEGRAL
00183               for (int ii=i+1; ii < dimE; ii++) result[ind] += (*ipart)->GetDensity(indspat*dimE+ii) * gasdensity * spall_spectrum[ii];
00184             }
00185           }
00186         }
00187       }
00188     }
00189   } // else
00190 
00191   return result;
00192 }

void TParticle::Evolve ( vector< TParticle * >  part,
vector< TCREvolutorBasis * >  crev_,
TSpallationNetwork *  spnet 
)

Propagate nuclear density from sources.

Parameters:
part Array of already propagated nuclei, to obtain the secondary source contribution.
crev_ Array of solution algorithms, to use them in cascade.

Definition at line 103 of file particle.cc.

References _fInXSec, A, ComputeSecondarySource(), daughter, density, eloss, TEnergyLoss::GetDpdt(), TSpectrum::GetSpectrum(), TInelasticCrossSection::GetXSec(), issec, lifetime, sp, and Z.

00103                                                                                                           {
00104 
00105   TEnergyLoss el(*(eloss.front()));
00106   for (vector<TEnergyLoss*>::iterator i = eloss.begin()+1; i != eloss.end(); ++i) el += (*(*i));
00107 
00108   const vector<double> secsource = ComputeSecondarySource(part, spnet);
00109 
00110   for (vector<TCREvolutorBasis*>::iterator crev = crev_.begin(); crev != crev_.end(); ++crev) (*crev)->Run(density, _fInXSec->GetXSec(), el.GetDpdt(), secsource, sp->GetSpectrum(), double(A), double(Z), lifetime, daughter, bool(issec));
00111 
00112   return ;
00113 }

double TParticle::FindNormalization (  ) 

Find normalization factor, by comparing with observed proton flux.

Returns:
the normalization factor.
Warning:
Use GSL cspline interpolator.

Definition at line 72 of file particle.cc.

References _fGalaxy, C, density, dimr, dimz, Galaxy::GetCoordinates(), TGrid::GetEk(), TGrid::GetR(), TGrid::GetZ(), index(), kpc, Rmax, Rmin, robs, sp_ref_rig_el, sp_ref_rig_norm, spect_norm, spect_norm_el, uid, and zobs.

Referenced by DRAGON::Run().

00072                                     {
00073   
00074   if (!(uid == 1001 || uid == -1000)) {
00075     cerr << "Asking for normalization but no protons nor primary electrons found" << endl;
00076     return -1;
00077   }
00078 
00079   vector<double> r = _fGalaxy->GetCoordinates()->GetR();
00080   unsigned int irsun = (unsigned int) ((robs-Rmin)/(Rmax-Rmin)*(double)(dimr-1));   
00081 
00082   vector<double> z = _fGalaxy->GetCoordinates()->GetZ();
00083   vector<double> Ek = _fGalaxy->GetCoordinates()->GetEk();
00084   
00085   unsigned int izsun = (unsigned int) ((zobs-z.front())/(z.back()-z.front())*(double)(dimz-1));
00086   double r1 = (r[irsun+1]-robs)/(r[irsun+1]-r[irsun]);
00087   double r2 = (robs-r[irsun])/(r[irsun+1]-r[irsun]);
00088 
00089   int ilow = (uid == 1001) ? int(log(sp_ref_rig_norm/Ek[0])/log(Ek[1]/Ek[0])) : int(log(sp_ref_rig_el/Ek[0])/log(Ek[1]/Ek[0]));
00090 
00091   double valuelow = log10( C/4.0/M_PI*pow(kpc,-3.0)*1e4*(density[index(irsun,izsun,ilow)]*r1+density[index(irsun+1,izsun,ilow)]*r2) );
00092   double valuehigh = log10( C/4.0/M_PI*pow(kpc,-3.0)*1e4*(density[index(irsun,izsun,ilow+1)]*r1+density[index(irsun+1,izsun,ilow+1)]*r2) );
00093 
00094   double e1 = (uid==1001) ? log(Ek[ilow+1]/sp_ref_rig_norm)/log(Ek[ilow+1]/Ek[ilow]) : log(Ek[ilow+1]/sp_ref_rig_el)/log(Ek[ilow+1]/Ek[ilow]);
00095   double e2 = (uid==1001) ? log(sp_ref_rig_norm/Ek[ilow])/log(Ek[ilow+1]/Ek[ilow]) : log(sp_ref_rig_el/Ek[ilow])/log(Ek[ilow+1]/Ek[ilow]);
00096   
00097   double value = pow(10,valuelow*e1 + valuehigh*e2);
00098   double factor = (uid == 1001) ? spect_norm/value : spect_norm_el/value;
00099 
00100   return factor;
00101 }

int TParticle::GetA (  )  [inline]

Get Mass number.

Definition at line 105 of file particle.h.

References A.

Referenced by DRAGON::particle_modulated_spectrum().

const int TParticle::GetDaughter (  )  const [inline]

Get UID of daughter nucleus.

Definition at line 104 of file particle.h.

References daughter.

double TParticle::GetDensity ( int  i  )  [inline]

Return density at given position and energy in linear representation.

Parameters:
i Linearized index

Definition at line 100 of file particle.h.

References density.

vector<double>& TParticle::GetDensity (  )  [inline]

Get density/spectrum array.

Definition at line 99 of file particle.h.

References density.

double TParticle::GetDensity ( int  ir,
int  iz,
int  ip 
) [inline]

Return density at given position and energy. Wrap from matrix to linear representation.

Parameters:
ir Radial position
iz Vertical position
ip Energy position

Definition at line 97 of file particle.h.

References density, and index().

Referenced by DRAGON::particle_modulated_spectrum().

int TParticle::GetIsSec (  )  [inline]

Definition at line 107 of file particle.h.

References issec.

00107 { return issec; }

const double TParticle::GetLifetime (  )  const [inline]

Get Life time

Definition at line 103 of file particle.h.

References lifetime.

TSpectrum* TParticle::GetSpectrum (  )  [inline]

Get injection spectrum.

Definition at line 96 of file particle.h.

References sp.

int TParticle::GetUid (  )  [inline]

Get unique nucleus ID = 1000*Z+A

Definition at line 102 of file particle.h.

References uid.

int TParticle::GetZ (  )  [inline]

Get charge.

Definition at line 106 of file particle.h.

References Z.

Referenced by DRAGON::particle_modulated_spectrum().

int TParticle::index ( int  ir,
int  iz,
int  ip 
) [inline, protected]

Convert matrix to linear representation.

Returns:
(ir*dimz+iz)*dimE+ip
Parameters:
ir Radial position
iz Vertical position
ip Energy position

Definition at line 132 of file particle.h.

References dimE, and dimz.

Referenced by FindNormalization(), GetDensity(), Print(), and PrintSpectrum().

00132 { return (ir*dimz+iz)*dimE+ip; }

TParticle& TParticle::operator+= ( const TParticle &  addendum  )  [inline]

To add nuclei of the same species but with different origin (e.g. primary and secondary protons, to obtain normalization).

Definition at line 74 of file particle.h.

References density, and uid.

void TParticle::Print ( fitsfile *  output_ptr,
double  norm 
)

Print propagated density and spectrum to file.

< Print all the information relevant to that nucleus: charge, mass, source abundance, injection spectrum and propagated density.

Parameters:
output_ptr Output file
norm Normalization factor

Definition at line 194 of file particle.cc.

References _fGalaxy, A, C, density, dimE, dimr, dimz, Galaxy::GetInjSpectrum(), Galaxy::GetSourceAbundance(), index(), issec, kpc, uid, and Z.

00194                                                        { 
00195   int status = 0;
00196 
00197   const long naxis = 3;
00198   long size_axes[naxis] = {dimE,dimr,dimz};
00199   long nelements = size_axes[0]*size_axes[1]*size_axes[2];
00200   long fpixel = 1;
00201   int bitpix = FLOAT_IMG;
00202   if (fits_create_img(output_ptr, bitpix, naxis, size_axes, &status)) fits_report_error(stderr, status);
00203 
00204   double sab = _fGalaxy->GetSourceAbundance(uid);
00205   double injind = _fGalaxy->GetInjSpectrum(uid);
00206   if (fits_write_key(output_ptr, TINT,    "Z_",      &Z,                 NULL, &status)) fits_report_error(stderr, status);
00207   if (fits_write_key(output_ptr, TINT,    "A",       &A,                 NULL, &status)) fits_report_error(stderr, status);
00208   if (fits_write_key(output_ptr, TINT,    "Sec",     &issec,             NULL, &status)) fits_report_error(stderr, status);
00209   if (fits_write_key(output_ptr, TDOUBLE, "S_Ab",    &sab,               NULL, &status)) fits_report_error(stderr, status);
00210   if (fits_write_key(output_ptr, TDOUBLE, "SpecInd", &injind,            NULL, &status)) fits_report_error(stderr, status);
00211 
00212   int counter = 0;
00213   float* array = new float[nelements]; 
00214   double weight = 1;
00215   if (A != 0) weight = double(A);
00216   for (int l = 0; l < dimz; ++l) {
00217     for (int k = 0; k < dimr; ++k) {
00218       for (int j = 0; j < dimE; ++j) {
00219         array[counter] = float(weight*density[index(k,l,j)]*norm*C/4.0/M_PI*pow(kpc,-3.0)*1e4);
00220         counter++;
00221       }
00222     }
00223   }
00224 
00225   if (fits_write_img(output_ptr, TFLOAT, fpixel, nelements, array, &status)) fits_report_error(stderr, status);
00226 
00227   delete [] array;
00228 
00229   return ;
00230 }

void TParticle::PrintSpectrum ( fitsfile *  output_ptr,
double  norm 
)

Print propagated spectrum to file.

< Print all the information relevant to that nucleus: charge, mass, source abundance, injection spectrum and propagated spectrum at Sun position.

Parameters:
output_ptr Output file
norm Normalization factor

Definition at line 232 of file particle.cc.

References _fGalaxy, A, C, density, dimE, dimr, dimz, Galaxy::GetCoordinates(), Galaxy::GetInjSpectrum(), TGrid::GetR(), Galaxy::GetSourceAbundance(), TGrid::GetZ(), index(), issec, kpc, Rmax, Rmin, robs, uid, Z, and zobs.

00232                                                                { 
00233   int status = 0;
00234   const long naxis = 1;
00235   long size_axes[naxis] = {dimE};
00236   long nelements = size_axes[0];
00237   long fpixel = 1;
00238   int bitpix = FLOAT_IMG;
00239   if (fits_create_img(output_ptr, bitpix, naxis, size_axes, &status)) fits_report_error(stderr, status);
00240   
00241   double sab = (!issec) ? _fGalaxy->GetSourceAbundance(uid) : 0.0;
00242   double injind = _fGalaxy->GetInjSpectrum(uid);
00243   if (fits_write_key(output_ptr, TINT,    "Z_",      &Z,                 NULL, &status)) fits_report_error(stderr, status);
00244   if (fits_write_key(output_ptr, TINT,    "A",       &A,                 NULL, &status)) fits_report_error(stderr, status);
00245   if (fits_write_key(output_ptr, TINT,    "Sec",     &issec,             NULL, &status)) fits_report_error(stderr, status);
00246   if (fits_write_key(output_ptr, TDOUBLE, "S_Ab",    &sab,               NULL, &status)) fits_report_error(stderr, status);
00247   if (fits_write_key(output_ptr, TDOUBLE, "SpecInd", &injind,            NULL, &status)) fits_report_error(stderr, status);
00248 
00249   float* array = new float[nelements]; 
00250   double weight = 1;
00251 
00252   vector<double> r = _fGalaxy->GetCoordinates()->GetR();
00253   unsigned int irsun = (unsigned int) ((robs-Rmin)/(Rmax-Rmin)*(double)(dimr-1));   
00254 
00255   unsigned int izsun = (unsigned int) ((zobs-_fGalaxy->GetCoordinates()->GetZ().front())/(2.0*_fGalaxy->GetCoordinates()->GetZ().back())*(double)(dimz-1));
00256 
00257   double r1 = (r[irsun+1]-robs)/(r[irsun+1]-r[irsun]);
00258   double r2 = (robs-r[irsun])/(r[irsun+1]-r[irsun]);
00259   if (A != 0) weight = double(A);
00260   for (int j = 0; j < dimE; ++j) {
00261     array[j] = float(weight*norm*( density[index(irsun,izsun,j)]*r1 + density[index(irsun+1,izsun,j)]*r2 )*C/4.0/M_PI*pow(kpc,-3.0)*1e4);
00262   }
00263   
00264   if (fits_write_img(output_ptr, TFLOAT, fpixel, nelements, array, &status)) fits_report_error(stderr, status);
00265   delete [] array;
00266   return ;
00267 }

Member Data Documentation

Galaxy* TParticle::_fGalaxy [protected]

galaxy model

Definition at line 119 of file particle.h.

Referenced by ComputeSecondarySource(), FindNormalization(), Print(), PrintSpectrum(), and TParticle().

TInelasticCrossSection* TParticle::_fInXSec [protected]

Inelastic cross section

Definition at line 120 of file particle.h.

Referenced by Evolve(), TParticle(), and ~TParticle().

int TParticle::A [protected]

Mass number

Definition at line 110 of file particle.h.

Referenced by ComputeSecondarySource(), Evolve(), GetA(), Print(), PrintSpectrum(), and TParticle().

int TParticle::daughter [protected]

UID of daughter nucleus

Definition at line 117 of file particle.h.

Referenced by Evolve(), GetDaughter(), and TParticle().

vector<double> TParticle::density [protected]

Spatial-energy distribution of propagated nucleus.

Definition at line 123 of file particle.h.

Referenced by ComputeSecondarySource(), Evolve(), FindNormalization(), GetDensity(), operator+=(), Print(), PrintSpectrum(), TParticle(), and ~TParticle().

int TParticle::dimE [protected]

energy dimension of the box

Definition at line 115 of file particle.h.

Referenced by ComputeSecondarySource(), index(), Print(), PrintSpectrum(), and TParticle().

int TParticle::dimr [protected]

radial dimension of the box

Definition at line 113 of file particle.h.

Referenced by ComputeSecondarySource(), FindNormalization(), Print(), PrintSpectrum(), and TParticle().

int TParticle::dimz [protected]

vertical dimension of the box

Definition at line 114 of file particle.h.

Referenced by ComputeSecondarySource(), FindNormalization(), index(), Print(), PrintSpectrum(), and TParticle().

vector<TEnergyLoss*> TParticle::eloss [protected]

Energy losses (several components)

Definition at line 122 of file particle.h.

Referenced by Evolve(), TParticle(), and ~TParticle().

int TParticle::issec [protected]

Whether is a secondary species (already exist a primary one).

Definition at line 118 of file particle.h.

Referenced by ComputeSecondarySource(), Evolve(), GetIsSec(), Print(), PrintSpectrum(), and TParticle().

double TParticle::lifetime [protected]

life time for unstable nuclei (=-1 if stable)

Definition at line 116 of file particle.h.

Referenced by Evolve(), GetLifetime(), and TParticle().

TSpectrum* TParticle::sp [protected]

Injection spectrum

Definition at line 121 of file particle.h.

Referenced by Evolve(), GetSpectrum(), TParticle(), and ~TParticle().

int TParticle::uid [protected]

Unique ID: = 1000*Z+A

Definition at line 112 of file particle.h.

Referenced by ComputeSecondarySource(), FindNormalization(), GetUid(), operator+=(), Print(), PrintSpectrum(), and TParticle().

int TParticle::Z [protected]

Charge

Definition at line 111 of file particle.h.

Referenced by Evolve(), GetZ(), Print(), PrintSpectrum(), and TParticle().

The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Enumerations Enumerator
Generated on Mon Sep 27 12:59:58 2010 for DRAGON by  doxygen 1.6.3