DRAGON Class Reference

DRAGON main class. More...

#include <dragon.h>

List of all members.

Public Member Functions

 DRAGON ()
 DRAGON (Input *in_)
 ~DRAGON ()
void Run (double &, double &)
void Print ()
void PrintChi2 ()
void PrintRatios (double modpotential=0.0)
 Print B/C, N/O, C/O and 10Be/9Be, for MCMC runs.

Protected Member Functions

string make_filename (const char *argv, const char *mode)
 Produce file name from run name.
vector< double > particle_modulated_spectrum (TParticle *particle, double potential)
 Modulate particle spectra at Sun position.
TParticle * FindParticle (int uid)
 Find particle in particle array.

Protected Attributes

Galaxy * gal
Input * in
vector< TCREvolutorBasis * > alg
TNucleiList * list
vector< TParticle * > particles
fitsfile * output_ptr
fitsfile * output_ptr_sp
int status

Detailed Description

DRAGON main class.

Author:
Luca Maccione luca.maccione@desy.de This class takes care of reading user input and settings, creating the necessary galactic structure, the algorithm to solve the diffusion equation, the structure to contain the propagated nuclear densities and the output.
See also:
Input
constants.h
TNucleiList
Galaxy
TCREvolutorBasis
TParticle

Definition at line 37 of file dragon.h.

Constructor & Destructor Documentation

DRAGON::DRAGON (  )  [inline]

The default constructor

Definition at line 40 of file dragon.h.

DRAGON::DRAGON ( Input *  in_  ) 

Constructor given some input.

Parameters:
in_ Wrapper to user input.

Definition at line 20 of file dragon.cc.

References ADI, alg, Input::filename, fullstore, gal, in, list, make_filename(), OpSplit, output_ptr, output_ptr_sp, partialstore, and status.

00020                                                         {
00021         in = in_; 
00022         
00023         cout << "Initializing the Galaxy...";
00024         gal = new Galaxy(in); 
00025         cout << "done" << endl;
00026         
00027         cout << "Initializing Algorithm...";
00028 
00029         if(OpSplit) alg.push_back( new TCREvolutor(gal) ); 
00030         if (ADI) alg.push_back( new TCREvolutorADI(gal) ); 
00031 
00032         cout << "done" << endl;
00033         
00034         cout << "Getting nuclei list... ";
00035         list = new TNucleiList();
00036         cout << " done" << endl;  
00037         
00038         if (fullstore) {
00039                 string name = make_filename(in->filename.c_str(), ".fits.gz");
00040                 if (fits_create_file(&output_ptr, name.c_str(), &status)) fits_report_error(stderr, status);
00041         }
00042         else output_ptr = NULL;
00043         
00044         if (partialstore) {
00045                 string namesp = make_filename(in->filename.c_str(), "_spectrum.fits.gz");
00046                 if (fits_create_file(&output_ptr_sp, namesp.c_str(), &status)) fits_report_error(stderr, status);
00047         }
00048         else output_ptr_sp = NULL;
00049         
00050 }

DRAGON::~DRAGON (  ) 

Destructor, that takes care also of closing the output files.

Definition at line 52 of file dragon.cc.

References alg, gal, in, list, output_ptr, output_ptr_sp, particles, and status.

00052                 {
00053         
00054         if (gal) delete gal;
00055         if (in) delete in;
00056         if (list) delete list;
00057         
00058         for (vector<TCREvolutorBasis*>::iterator i = alg.begin(); i != alg.end(); ++i) delete *i;
00059         alg.clear();
00060         
00061         for (vector<TParticle*>::iterator i = particles.begin(); i != particles.end(); ++i) delete *i;
00062         particles.clear();
00063         
00064         if (output_ptr) {
00065                 if (fits_close_file(output_ptr, &status)) fits_report_error(stderr, status);
00066         }
00067         if (output_ptr_sp) {
00068                 if (fits_close_file(output_ptr_sp, &status)) fits_report_error(stderr, status);
00069         }
00070 }

Member Function Documentation

TParticle * DRAGON::FindParticle ( int  uid  )  [protected]

Find particle in particle array.

Parameters:
uid Unique ID of nuclear species.
Returns:
Particle density of specified type

Definition at line 72 of file dragon.cc.

References particles.

Referenced by PrintChi2(), and PrintRatios().

00072                                        {
00073         for (vector<TParticle*>::reverse_iterator ripart = particles.rbegin(); ripart != particles.rend(); ++ripart) {
00074                 if ( (*ripart)->GetUid() == uid ) return *ripart;
00075         }
00076         
00077         return NULL;
00078 }

string DRAGON::make_filename ( const char *  argv,
const char *  mode 
) [inline, protected]

Produce file name from run name.

Parameters:
argv Run name
mode extension

Definition at line 72 of file dragon.h.

Referenced by DRAGON().

00072                                                            {
00073     string filename = "!output/";
00074     filename += argv;   
00075     return filename + mode;
00076   }

vector< double > DRAGON::particle_modulated_spectrum ( TParticle *  particle,
double  potential 
) [protected]

Modulate particle spectra at Sun position.

Parameters:
particle Pointer to a nuclear species.
potential Modulation potential.
Returns:
Modulated spectrum at Sun position.
Warning:
Make use of GSL cspline.

Definition at line 80 of file dragon.cc.

References ADI, gal, TParticle::GetA(), Galaxy::GetCoordinates(), TGrid::GetDeltaR(), TGrid::GetDeltaZ(), TParticle::GetDensity(), TGrid::GetDimE(), TGrid::GetEk(), TGrid::GetR(), TParticle::GetZ(), TGrid::GetZ(), mp, robs, and zobs.

Referenced by PrintChi2(), and PrintRatios().

00080                                                                                         {
00081         
00082         if (particle == NULL) {
00083                 cerr << "Null Particle in particle_modulated_spectrum" << endl;
00084                 //exit(3);
00085                 return vector<double>(gal->GetCoordinates()->GetDimE(),0.0);
00086         }
00087         
00088         // Sun position
00089         vector<double> r = gal->GetCoordinates()->GetR();
00090         int irsun = (int) (robs/gal->GetCoordinates()->GetDeltaR()); 
00091 
00092         int izsun = (int) ((zobs-gal->GetCoordinates()->GetZ().front())/gal->GetCoordinates()->GetDeltaZ());
00093         
00094         int dimE = gal->GetCoordinates()->GetDimE();
00095         vector<double> Ekin = gal->GetCoordinates()->GetEk();
00096         
00097         double r1 = (r[irsun+1]-robs)/(r[irsun+1]-r[irsun]);
00098         double r2 = (robs-r[irsun])/(r[irsun+1]-r[irsun]);
00099         double spectrum = 0.;
00100         
00101         double logEk[dimE-(ADI)*2];
00102         double particlesp1[dimE-(ADI)*2];
00103         for (int m = 0+(ADI); m < dimE-(ADI); m++) {
00104           spectrum = (particle->GetDensity(irsun,izsun,m)*r1 + particle->GetDensity(irsun+1,izsun,m)*r2)*double(particle->GetA());
00105           particlesp1[m-(ADI)] = (spectrum != 0.) ? log10(spectrum) : 100.;
00106           logEk[m-(ADI)] = log10(Ekin[m]);
00107         }
00108         
00109         double Phi = double(particle->GetZ())/double(particle->GetA())*potential;
00110         double T0 = 0.;
00111         gsl_interp_accel *acc = gsl_interp_accel_alloc();
00112         gsl_spline *particlespline = gsl_spline_alloc(gsl_interp_cspline, dimE-(ADI)*2);
00113         gsl_spline_init(particlespline, logEk, particlesp1, dimE-(ADI)*2);
00114         
00115         vector<double> particlesp(dimE, 0.0);
00116         for (int m = 0+(ADI); m < dimE-(ADI); m++) {            
00117                 T0 = Ekin[m] + Phi;
00118                 particlesp[m] = pow(10.0, gsl_spline_eval(particlespline, log10(T0), acc))*(Ekin[m]*(Ekin[m]+2.0*mp))/(T0*(T0+2.0*mp));
00119         }
00120         
00121         gsl_spline_free(particlespline);
00122         gsl_interp_accel_free(acc);
00123         
00124         return particlesp;
00125 }

void DRAGON::Print (  ) 

Print the first HDU, with general information on the run.

< Radial Sun position.

< Vertical Sun position.

Create the first HDU, which only contains header data.

Definition at line 127 of file dragon.cc.

References Bh, Input::D0, D_ref_rig, Input::delta, Input::dmmode, Input::dmprof, Input::dvdz, Input::Ekfact, Ekmax, Ekmin, gas_model, He_abundance, in, km, kpc, Input::mx, Myr, Input::numr, Input::numz, output_ptr, output_ptr_sp, p, rB, Rmax, Rmin, robs, Input::sigmav, SNR_model, sp_ref_rig, sp_ref_rig_norm, spect_norm, status, Input::taudec, u, Input::v0, Input::vAlfven, Input::zmax, zobs, zr, and Input::zt.

Referenced by main().

00127                    {
00128   unsigned int cr_irsun = (unsigned int) (robs/Rmax*(double)(in->numr-1)); 
00130   unsigned int cr_izsun = (unsigned int) ((zobs+in->zmax)/(2.0*in->zmax)*(double)(in->numz-1)), 
00131     cr_gas_model = (unsigned int)gas_model,
00132     cr_SNR_model = (unsigned int)SNR_model;
00133   
00134   double cr_rmin = Rmin,
00135     cr_rmax = Rmax,
00136     cr_zmin = -in->zmax,
00137     cr_zmax = in->zmax,
00138     cr_rB = rB,
00139     cr_zt = in->zt,
00140     cr_zr = zr,
00141     cr_Bh = Bh,
00142     cr_D0 = in->D0*kpc*kpc/Myr,
00143     cr_D_ref_rig = D_ref_rig,
00144     cr_ind_diff = in->delta,
00145     cr_He_abundance = He_abundance,
00146     cr_sp_ref_rig = sp_ref_rig,
00147     cr_sp_ref_rig_norm = sp_ref_rig_norm,
00148     cr_spect_norm = spect_norm,
00149     cr_Ekmax = Ekmax,
00150     cr_Ekmin = Ekmin,
00151     cr_Ekfac = in->Ekfact,
00152     cr_u = u,
00153     cr_p = p,
00154     cr_robs = robs,
00155 #ifdef REAC
00156     cr_vA = in->vAlfven*kpc/Myr/km,
00157 #endif
00158 #ifdef CONVECTION
00159     cr_v0 = in->v0*kpc/Myr/km,
00160     cr_dvdz = in->dvdz*kpc/Myr/km,
00161 #endif
00162     cr_zobs = zobs;
00163 
00164 #ifdef HAVE_DS
00165   double crmx = in->mx;
00166   double crtaudec = in->taudec;
00167   double crsigmav = in->sigmav;
00168   int crdmmode = in->dmmode;
00169   int crdmprof = in->dmprof;
00170 #endif
00171 
00172   const long naxis = 3;
00173   int dimE = int(log(Ekmax/Ekmin)/log(in->Ekfact) + 1.9);
00174   long size_axes[naxis] = {dimE,in->numr,in->numz};
00175   long nelements = size_axes[0]*size_axes[1]*size_axes[2];
00176   long fpixel = 1;
00177   int bitpix = FLOAT_IMG;
00178 
00183   if (output_ptr) {
00184     if (fits_create_img(output_ptr, bitpix, naxis, size_axes, &status)) fits_report_error(stderr, status);
00185     
00186     if (fits_write_key(output_ptr, TDOUBLE, "rmin",      &cr_rmin,            NULL, &status)) fits_report_error(stderr, status);
00187     if (fits_write_key(output_ptr, TDOUBLE, "rmax",      &cr_rmax,            NULL, &status)) fits_report_error(stderr, status);
00188     if (fits_write_key(output_ptr, TDOUBLE, "zmin",      &cr_zmin,            NULL, &status)) fits_report_error(stderr, status);
00189     if (fits_write_key(output_ptr, TDOUBLE, "zmax",      &cr_zmax,            NULL, &status)) fits_report_error(stderr, status);
00190     if (fits_write_key(output_ptr, TDOUBLE, "robs",      &cr_robs,            NULL, &status)) fits_report_error(stderr, status); 
00191     if (fits_write_key(output_ptr, TDOUBLE, "zobs",      &cr_zobs,            NULL, &status)) fits_report_error(stderr, status);
00192     if (fits_write_key(output_ptr, TUINT,   "irsun",     &cr_irsun,           NULL, &status)) fits_report_error(stderr, status); 
00193     if (fits_write_key(output_ptr, TUINT,   "izsun",     &cr_izsun,           NULL, &status)) fits_report_error(stderr, status);
00194     if (fits_write_key(output_ptr, TDOUBLE, "u",         &cr_u,               NULL, &status)) fits_report_error(stderr, status);
00195     if (fits_write_key(output_ptr, TDOUBLE, "rB",        &cr_rB,              NULL, &status)) fits_report_error(stderr, status);
00196     if (fits_write_key(output_ptr, TDOUBLE, "zt",        &cr_zt,              NULL, &status)) fits_report_error(stderr, status);
00197     if (fits_write_key(output_ptr, TDOUBLE, "zr",        &cr_zr,              NULL, &status)) fits_report_error(stderr, status);
00198     if (fits_write_key(output_ptr, TDOUBLE, "Bh",        &cr_Bh,              NULL, &status)) fits_report_error(stderr, status);
00199     if (fits_write_key(output_ptr, TDOUBLE, "D0",        &cr_D0,              NULL, &status)) fits_report_error(stderr, status);
00200     if (fits_write_key(output_ptr, TDOUBLE, "Drefrig",   &cr_D_ref_rig,       NULL, &status)) fits_report_error(stderr, status);
00201     if (fits_write_key(output_ptr, TDOUBLE, "ind_diff",  &cr_ind_diff,        NULL, &status)) fits_report_error(stderr, status);
00202     if (fits_write_key(output_ptr, TDOUBLE, "He_ab",     &cr_He_abundance,    NULL, &status)) fits_report_error(stderr, status);
00203     if (fits_write_key(output_ptr, TDOUBLE, "sprefrig",  &cr_sp_ref_rig,      NULL, &status)) fits_report_error(stderr, status);
00204     if (fits_write_key(output_ptr, TDOUBLE, "sprignor",  &cr_sp_ref_rig_norm, NULL, &status)) fits_report_error(stderr, status);
00205     if (fits_write_key(output_ptr, TDOUBLE, "specnorm",  &cr_spect_norm,      NULL, &status)) fits_report_error(stderr, status);
00206     if (fits_write_key(output_ptr, TDOUBLE, "Ekmax",     &cr_Ekmax,           NULL, &status)) fits_report_error(stderr, status);
00207     if (fits_write_key(output_ptr, TDOUBLE, "Ekmin",     &cr_Ekmin,           NULL, &status)) fits_report_error(stderr, status);
00208     if (fits_write_key(output_ptr, TDOUBLE, "Ekin_fac",  &cr_Ekfac,           NULL, &status)) fits_report_error(stderr, status);
00209     if (fits_write_key(output_ptr, TUINT,   "gasmod",    &cr_gas_model,       NULL, &status)) fits_report_error(stderr, status);
00210     if (fits_write_key(output_ptr, TUINT,   "SNRmod",    &cr_SNR_model,       NULL, &status)) fits_report_error(stderr, status);
00211     if (fits_write_key(output_ptr, TINT,    "dimE",      &dimE,               NULL, &status)) fits_report_error(stderr, status);
00212     if (fits_write_key(output_ptr, TINT,    "dimr",      &(in->numr),         NULL, &status)) fits_report_error(stderr, status);
00213     if (fits_write_key(output_ptr, TINT,    "dimz",      &(in->numz),         NULL, &status)) fits_report_error(stderr, status);
00214 #ifdef REAC
00215     if (fits_write_key(output_ptr, TDOUBLE, "vAlfven",   &cr_vA,              NULL, &status)) fits_report_error(stderr, status);
00216 #endif        
00217 #ifdef CONVECTION
00218     if (fits_write_key(output_ptr, TDOUBLE, "v0_conv",   &cr_v0,              NULL, &status)) fits_report_error(stderr, status);
00219     if (fits_write_key(output_ptr, TDOUBLE, "dvdz_conv", &cr_dvdz,            NULL, &status)) fits_report_error(stderr, status);
00220 #endif
00221 #ifdef HAVE_DS
00222     if (fits_write_key(output_ptr, TDOUBLE, "mx",        &crmx,               NULL, &status)) fits_report_error(stderr, status);
00223     if (fits_write_key(output_ptr, TDOUBLE, "taudec",    &crtaudec,           NULL, &status)) fits_report_error(stderr, status);
00224     if (fits_write_key(output_ptr, TDOUBLE, "sigmav",    &crsigmav,           NULL, &status)) fits_report_error(stderr, status);
00225     if (fits_write_key(output_ptr, TINT,    "dmmode",    &crdmmode,           NULL, &status)) fits_report_error(stderr, status);
00226     if (fits_write_key(output_ptr, TINT,    "dmprof",    &crdmprof,           NULL, &status)) fits_report_error(stderr, status);
00227 #endif
00228   }
00229   if (output_ptr_sp) {
00230     if (fits_create_img(output_ptr_sp, bitpix, naxis, size_axes, &status)) fits_report_error(stderr, status);
00231     
00232     if (fits_write_key(output_ptr_sp, TDOUBLE, "rmin",      &cr_rmin,            NULL, &status)) fits_report_error(stderr, status);
00233     if (fits_write_key(output_ptr_sp, TDOUBLE, "rmax",      &cr_rmax,            NULL, &status)) fits_report_error(stderr, status);
00234     if (fits_write_key(output_ptr_sp, TDOUBLE, "zmin",      &cr_zmin,            NULL, &status)) fits_report_error(stderr, status);
00235     if (fits_write_key(output_ptr_sp, TDOUBLE, "zmax",      &cr_zmax,            NULL, &status)) fits_report_error(stderr, status);
00236     if (fits_write_key(output_ptr_sp, TDOUBLE, "robs",      &cr_robs,            NULL, &status)) fits_report_error(stderr, status); 
00237     if (fits_write_key(output_ptr_sp, TDOUBLE, "zobs",      &cr_zobs,            NULL, &status)) fits_report_error(stderr, status);
00238     if (fits_write_key(output_ptr_sp, TUINT,   "irsun",     &cr_irsun,           NULL, &status)) fits_report_error(stderr, status); 
00239     if (fits_write_key(output_ptr_sp, TUINT,   "izsun",     &cr_izsun,           NULL, &status)) fits_report_error(stderr, status);
00240     if (fits_write_key(output_ptr_sp, TDOUBLE, "u",         &cr_u,               NULL, &status)) fits_report_error(stderr, status);
00241     if (fits_write_key(output_ptr_sp, TDOUBLE, "rB",        &cr_rB,              NULL, &status)) fits_report_error(stderr, status);
00242     if (fits_write_key(output_ptr_sp, TDOUBLE, "zt",        &cr_zt,              NULL, &status)) fits_report_error(stderr, status);
00243     if (fits_write_key(output_ptr_sp, TDOUBLE, "zr",        &cr_zr,              NULL, &status)) fits_report_error(stderr, status);
00244     if (fits_write_key(output_ptr_sp, TDOUBLE, "Bh",        &cr_Bh,              NULL, &status)) fits_report_error(stderr, status);
00245     if (fits_write_key(output_ptr_sp, TDOUBLE, "D0",        &cr_D0,              NULL, &status)) fits_report_error(stderr, status);
00246     if (fits_write_key(output_ptr_sp, TDOUBLE, "Drefrig",   &cr_D_ref_rig,       NULL, &status)) fits_report_error(stderr, status);
00247     if (fits_write_key(output_ptr_sp, TDOUBLE, "ind_diff",  &cr_ind_diff,        NULL, &status)) fits_report_error(stderr, status);
00248     if (fits_write_key(output_ptr_sp, TDOUBLE, "He_ab",     &cr_He_abundance,    NULL, &status)) fits_report_error(stderr, status);
00249     if (fits_write_key(output_ptr_sp, TDOUBLE, "sprefrig",  &cr_sp_ref_rig,      NULL, &status)) fits_report_error(stderr, status);
00250     if (fits_write_key(output_ptr_sp, TDOUBLE, "sprignor",  &cr_sp_ref_rig_norm, NULL, &status)) fits_report_error(stderr, status);
00251     if (fits_write_key(output_ptr_sp, TDOUBLE, "specnorm",  &cr_spect_norm,      NULL, &status)) fits_report_error(stderr, status);
00252     if (fits_write_key(output_ptr_sp, TDOUBLE, "Ekmax",     &cr_Ekmax,           NULL, &status)) fits_report_error(stderr, status);
00253     if (fits_write_key(output_ptr_sp, TDOUBLE, "Ekmin",     &cr_Ekmin,           NULL, &status)) fits_report_error(stderr, status);
00254     if (fits_write_key(output_ptr_sp, TDOUBLE, "Ekin_fac",  &cr_Ekfac,           NULL, &status)) fits_report_error(stderr, status);
00255     //   if (fits_write_key(output_ptr_sp, TDOUBLE, "tol",       &cr_tol,             NULL, &status)) fits_report_error(stderr, status);
00256     if (fits_write_key(output_ptr_sp, TUINT,   "gasmod",    &cr_gas_model,       NULL, &status)) fits_report_error(stderr, status);
00257     if (fits_write_key(output_ptr_sp, TUINT,   "SNRmod",    &cr_SNR_model,       NULL, &status)) fits_report_error(stderr, status);
00258     if (fits_write_key(output_ptr_sp, TINT,    "dimE",      &dimE,               NULL, &status)) fits_report_error(stderr, status);
00259     if (fits_write_key(output_ptr_sp, TINT,    "dimr",      &(in->numr),         NULL, &status)) fits_report_error(stderr, status);
00260     if (fits_write_key(output_ptr_sp, TINT,    "dimz",      &(in->numz),         NULL, &status)) fits_report_error(stderr, status);
00261 #ifdef REAC
00262     if (fits_write_key(output_ptr_sp, TDOUBLE, "vAlfven",   &cr_vA,              NULL, &status)) fits_report_error(stderr, status);
00263 #endif        
00264 #ifdef CONVECTION
00265     if (fits_write_key(output_ptr_sp, TDOUBLE, "v0_conv",   &cr_v0,              NULL, &status)) fits_report_error(stderr, status);
00266     if (fits_write_key(output_ptr_sp, TDOUBLE, "dvdz_conv", &cr_dvdz,            NULL, &status)) fits_report_error(stderr, status);
00267 #endif
00268 #ifdef HAVE_DS
00269     if (fits_write_key(output_ptr_sp, TDOUBLE, "mx",        &crmx,               NULL, &status)) fits_report_error(stderr, status);
00270     if (fits_write_key(output_ptr_sp, TDOUBLE, "taudec",    &crtaudec,           NULL, &status)) fits_report_error(stderr, status);
00271     if (fits_write_key(output_ptr_sp, TDOUBLE, "sigmav",    &crsigmav,           NULL, &status)) fits_report_error(stderr, status);
00272     if (fits_write_key(output_ptr_sp, TINT,    "dmmode",    &crdmmode,           NULL, &status)) fits_report_error(stderr, status);
00273     if (fits_write_key(output_ptr_sp, TINT,    "dmprof",    &crdmprof,           NULL, &status)) fits_report_error(stderr, status);
00274 #endif
00275   }
00276     
00277   float* array = new float[nelements]();
00278   if (output_ptr) {
00279     if (fits_write_img(output_ptr, TFLOAT, fpixel, nelements, array, &status)) fits_report_error(stderr, status);
00280   }
00281   if (output_ptr_sp) {
00282     if (fits_write_img(output_ptr_sp, TFLOAT, fpixel, nelements, array, &status)) fits_report_error(stderr, status);
00283   }
00284   delete [] array;
00285 }

void DRAGON::PrintChi2 (  ) 

Compute and print chi2, for MCMC runs

Definition at line 376 of file dragon.cc.

References Utility::chi2_struct::all, Utility::chi2_struct::ATIC, Utility::chi2_struct::CREAM, Utility::dimDATA, Utility::dimDATABe, Input::filename, FindParticle(), gal, Galaxy::GetCoordinates(), TGrid::GetDimE(), TGrid::GetEk(), Utility::chi2_struct::HEAO3, in, Utility::chi2_struct::N_all, Utility::chi2_struct::N_ATIC, Utility::chi2_struct::N_CREAM, Utility::chi2_struct::N_HEAO3, and particle_modulated_spectrum().

00376                        {
00377   
00378   data data_setBC[dimDATA];
00379   data_be data_setBe[dimDATABe];
00380 
00381 #include "chi2data.h"
00382 
00383   string name = "output/";
00384   name += in->filename.c_str();
00385   name += "_chi2mod.dat";
00386   
00387   ofstream outf;
00388   outf.open(name.c_str(), ios::out);
00389 
00390   int dimE = gal->GetCoordinates()->GetDimE();
00391   vector<double> Ekin(gal->GetCoordinates()->GetEk());
00392   
00393   double interpolated_Be = 0.0, 
00394     interpolated_BC = 0.0, 
00395     interpolated_CO = 0.0, 
00396     interpolated_NO = 0.0,
00397     delta_Be = 0.0,
00398     delta_BC = 0.0, 
00399     delta_CO = 0.0, 
00400     delta_NO = 0.0;
00401   
00402   chi2_struct chi2 = {0.0,0.0,0.0,0,
00403                       0.0,0.0,0.0,0,
00404                       0.0,0.0,0.0,0,
00405                       0.0,0.0,0.0,0};
00406   
00407   double chi2Be = 0.0;
00408   int chi2_NBe  = 0;
00409   
00410   double logBe[dimE]; 
00411   double logBC[dimE];
00412   double logCO[dimE];
00413   double logNO[dimE];
00414   double logEk[dimE];
00415  
00416   for (int i=0; i<dimE; ++i)
00417     logEk[i] = log10(Ekin[i]);
00418  
00419   // CHI2 WITH BE
00420   for (int j = 0; j < dimDATABe; ++j) {
00421     if (data_setBe[j].E > Ekin[0] && data_setBe[j].E < Ekin[dimE-1]) {
00422       
00423       vector<double> Be9  = particle_modulated_spectrum(FindParticle(4009), data_setBe[j].phi);
00424       vector<double> Be10 = particle_modulated_spectrum(FindParticle(4010), data_setBe[j].phi);
00425       
00426       for (int i = 0; i < dimE; ++i)
00427         logBe[i] = (Be9[i]!=0.) ? log10(Be10[i]/Be9[i]) : 100.0;
00428       
00429       gsl_interp_accel *acc = gsl_interp_accel_alloc ();
00430       gsl_spline *bespline = gsl_spline_alloc (gsl_interp_cspline, dimE);
00431       gsl_spline_init (bespline, logEk, logBe, dimE);
00432       
00433       interpolated_Be = gsl_spline_eval (bespline, log10(data_setBe[j].E), acc);
00434       
00435       delta_Be  = pow(data_setBe[j].Be - pow(10.0, interpolated_Be),2.0);
00436       delta_Be /= pow(data_setBe[j].sigma_Be,2.0);
00437       
00438       chi2Be += delta_Be;
00439       chi2_NBe++;
00440       
00441       gsl_spline_free (bespline);
00442       gsl_interp_accel_free (acc);
00443     }
00444   }
00445 
00446   // CHI2 WITH NUCLEI
00447   for (int j = 0; j < dimDATA; ++j) {
00448     if (data_setBC[j].E > Ekin[0] && data_setBC[j].E < Ekin[dimE-1]) {
00449 
00450       vector<double> B10 = particle_modulated_spectrum(FindParticle(5010), data_setBC[j].phi);
00451       vector<double> B11 = particle_modulated_spectrum(FindParticle(5011), data_setBC[j].phi);
00452       
00453       vector<double> C12 = particle_modulated_spectrum(FindParticle(6012), data_setBC[j].phi);
00454       vector<double> C13 = particle_modulated_spectrum(FindParticle(6013), data_setBC[j].phi);
00455       vector<double> C14 = particle_modulated_spectrum(FindParticle(6014), data_setBC[j].phi);
00456       
00457       vector<double> N14 = particle_modulated_spectrum(FindParticle(7014), data_setBC[j].phi);
00458       vector<double> N15 = particle_modulated_spectrum(FindParticle(7015), data_setBC[j].phi); 
00459       
00460       vector<double> O16 = particle_modulated_spectrum(FindParticle(8016), data_setBC[j].phi);
00461       vector<double> O17 = particle_modulated_spectrum(FindParticle(8017), data_setBC[j].phi);
00462       vector<double> O18 = particle_modulated_spectrum(FindParticle(8018), data_setBC[j].phi);
00463       
00464       for (int i = 0; i < dimE; ++i) {
00465         logBC[i] = (C12[i]+C13[i]+C14[i] !=0.) ? log10((B10[i]+B11[i])/(C12[i]+C13[i]+C14[i])) : 100.0;
00466         logCO[i] = (O16[i]+O17[i]+O18[i] !=0.) ? log10((C12[i]+C13[i])/(O16[i]+O17[i]+O18[i])) : 100.0;
00467         logNO[i] = (O16[i]+O17[i]+O18[i] !=0.) ? log10((N14[i]+N15[i])/(O16[i]+O17[i]+O18[i])) : 100.0;
00468         logEk[i] = log10(Ekin[i]);
00469       }
00470      
00471       gsl_interp_accel *acc = gsl_interp_accel_alloc ();
00472       gsl_spline *bcspline = gsl_spline_alloc (gsl_interp_cspline, dimE);
00473       gsl_spline *cospline = gsl_spline_alloc (gsl_interp_cspline, dimE);
00474       gsl_spline *nospline = gsl_spline_alloc (gsl_interp_cspline, dimE);
00475       
00476       gsl_spline_init (bcspline, logEk, logBC, dimE);
00477       gsl_spline_init (cospline, logEk, logCO, dimE);
00478       gsl_spline_init (nospline, logEk, logNO, dimE);
00479       
00480       interpolated_BC = gsl_spline_eval (bcspline, log10(data_setBC[j].E), acc);
00481       interpolated_CO = gsl_spline_eval (cospline, log10(data_setBC[j].E), acc);
00482       interpolated_NO = gsl_spline_eval (nospline, log10(data_setBC[j].E), acc);
00483       
00484       delta_BC  = pow(data_setBC[j].BC - pow(10.0, interpolated_BC),2.0);
00485       delta_BC /= pow(data_setBC[j].sigma_BC,2.0);
00486       delta_CO  = pow(data_setBC[j].CO - pow(10.0, interpolated_CO),2.0);
00487       delta_CO /= pow(data_setBC[j].sigma_CO,2.0);
00488       delta_NO  = pow(data_setBC[j].NO - pow(10.0, interpolated_NO),2.0);
00489       delta_NO /= pow(data_setBC[j].sigma_NO,2.0);
00490       
00491       chi2.all[0] += delta_BC;
00492       chi2.all[1] += delta_CO;
00493       chi2.all[2] += delta_NO;
00494       chi2.N_all++;
00495       
00496       if (data_setBC[j].experiment == 0) {
00497         chi2.HEAO3[0] += delta_BC;
00498         chi2.HEAO3[1] += delta_CO;
00499         chi2.HEAO3[2] += delta_NO;
00500         chi2.N_HEAO3++;
00501       }
00502       if (data_setBC[j].experiment == 3) {
00503         chi2.CREAM[0] += delta_BC;
00504         chi2.CREAM[1] += delta_CO;
00505         chi2.CREAM[2] += delta_NO;
00506         chi2.N_CREAM++;
00507       }
00508       if (data_setBC[j].experiment == 4) {
00509         chi2.ATIC[0] += delta_BC;
00510         chi2.ATIC[1] += delta_CO;
00511         chi2.ATIC[2] += delta_NO;
00512         chi2.N_ATIC++;
00513       }
00514       
00515       gsl_spline_free (bcspline);
00516       gsl_spline_free (cospline);
00517       gsl_spline_free (nospline);
00518       gsl_interp_accel_free (acc);
00519     }
00520   }
00521   
00522   for (int i = 0; i < 3; i++) {
00523     chi2.all[i]   /= (double)chi2.N_all-3;
00524     chi2.HEAO3[i] /= (double)chi2.N_HEAO3-3;
00525     chi2.CREAM[i] /= (double)chi2.N_CREAM-3;
00526     chi2.ATIC[i]  /= (double)chi2.N_ATIC-3;
00527   }
00528   
00529   outf<<"#ALL\t"<<"HEAO3\t"<<"CREAM\t"<<"ATIC\t"<<endl;
00530   for (int i = 0; i < 3; i++) 
00531       outf << chi2.all[i] << "\t" << chi2.HEAO3[i] << "\t" << chi2.CREAM[i] << "\t" << chi2.ATIC[i] << "\t" << endl;
00532 
00533   outf << chi2Be/(double)chi2_NBe << endl;
00534   outf.close();
00535   return;
00536 }

void DRAGON::PrintRatios ( double  modpotential = 0.0  ) 

Print B/C, N/O, C/O and 10Be/9Be, for MCMC runs.

Parameters:
modpotential [GV] Modulation potential

Definition at line 538 of file dragon.cc.

References ADI, Input::filename, FindParticle(), gal, Galaxy::GetCoordinates(), TGrid::GetDimE(), TGrid::GetEk(), in, and particle_modulated_spectrum().

00538                                             {
00539 
00540   string name = "output/";
00541   name += in->filename.c_str();
00542   name += "_ratios.dat";
00543   cout << "Printing ratios" << endl;
00544   ofstream outf;
00545   outf.open(name.c_str(), ios::out);
00546   
00547   int dimE = gal->GetCoordinates()->GetDimE();
00548   vector<double> Ekin(gal->GetCoordinates()->GetEk());
00549 
00550   vector<double> Be9  = particle_modulated_spectrum(FindParticle(4009), modpotential);
00551   vector<double> Be10 = particle_modulated_spectrum(FindParticle(4010), modpotential);
00552 
00553   vector<double> B10 = particle_modulated_spectrum(FindParticle(5010), modpotential);
00554   vector<double> B11 = particle_modulated_spectrum(FindParticle(5011), modpotential);
00555   
00556   vector<double> C12 = particle_modulated_spectrum(FindParticle(6012), modpotential);
00557   vector<double> C13 = particle_modulated_spectrum(FindParticle(6013), modpotential);
00558   vector<double> C14 = particle_modulated_spectrum(FindParticle(6014), modpotential);
00559   
00560   vector<double> N14 = particle_modulated_spectrum(FindParticle(7014), modpotential);
00561   vector<double> N15 = particle_modulated_spectrum(FindParticle(7015), modpotential); 
00562   
00563   vector<double> O16 = particle_modulated_spectrum(FindParticle(8016), modpotential);
00564   vector<double> O17 = particle_modulated_spectrum(FindParticle(8017), modpotential);
00565   vector<double> O18 = particle_modulated_spectrum(FindParticle(8018), modpotential);
00566 
00567   outf << "#Ek" << "\t" << "B/C" << "\t" << "C/O" << "\t" << "N/O" << "\t" << "Be10/Be9" << endl;
00568   for(int i = 0+(ADI); i < dimE-(ADI); ++i) {
00569 
00570     outf << Ekin[i] << "\t" << (B10[i]+B11[i])/(C12[i]+C13[i]+C14[i]) << "\t" << (C12[i]+C13[i]+C14[i])/(O16[i]+O17[i]+O18[i]) << "\t" << (N14[i]+N15[i])/(O16[i]+O17[i]+O18[i]) << "\t" << Be10[i]/Be9[i] << endl;
00571   }
00572   outf.close();
00573   return;
00574 }

void DRAGON::Run ( double &  norm,
double &  normel 
)

Main routine of DRAGON, that actually propagates the nuclear chain.

< Compute the ID of the nucleus.

< If nucleus has NO source abundance then it is not propagated. Just a way to control propagation of specific isotopes

< Creates the nucleus and add it to the particle vector.

< Evolve the nucleus.

Normalization part. If it is DM then just convert units, otherwise find protons and electrons and normalize to observed flux. Units are GeV^-1 m^-2 s^-1 sr^-1.

Print nuclear densities and/or spectra in output file. Lighter nuclei are printed first.

Definition at line 287 of file dragon.cc.

References alg, C, TParticle::FindNormalization(), gal, GalpropXSec, Galaxy::GetCoordinates(), TNucleiList::GetList(), Galaxy::GetSourceAbundance(), Utility::id_nuc(), in, kpc, list, output_ptr, output_ptr_sp, particles, Pi, prop_extracomp, spallationxsec, and Webber03.

Referenced by main().

00287                                              {
00288 
00289   norm = 1.0;
00290   normel = 1.0;
00291 
00292   TXSecBase* xsecmodel = NULL;
00293   if (spallationxsec == GalpropXSec) xsecmodel = new TGalpropXSec(gal->GetCoordinates());
00294   else if (spallationxsec == Webber03) xsecmodel = new TWebber03();
00295   else {
00296         cerr << "Problem with spallation xsec models. Exiting." << endl;
00297 #ifndef HAVE_DS
00298         exit(-4);
00299 #endif
00300   }
00301   vector<int> list_nuc = list->GetList();
00302 
00303   TSpallationNetwork* spallnet = new TSpallationNetwork(gal->GetCoordinates(), xsecmodel, list_nuc);
00304 
00305   for (vector<int>::iterator inuc = list_nuc.begin(); inuc != list_nuc.end(); ++inuc) {
00306 
00307     int A = -1000;
00308     int Z = -1000;
00309     Utility::id_nuc(*inuc, A, Z); 
00311     cout << "Starting with nucleus A = " << A << " Z = " << Z << endl;
00312     
00313     if (gal->GetSourceAbundance(*inuc) < 0) { 
00314       cerr << "Nucleus not included in initialization list. Not propagating...\n" << endl;
00315       continue;
00316     }
00317 
00318     cout << "Starting propagation..." << endl;
00319     int prev_uid = (particles.size() > 0) ? particles.back()->GetUid() : 0;
00320     particles.push_back(new TParticle(A, Z, gal, in, prev_uid, xsecmodel, list)); 
00321     particles.back()->Evolve(particles, alg, spallnet); 
00323     cout << "Propagation done.\n" << endl;
00324   }
00325 
00329 #ifdef HAVE_DS
00330   norm *= C * pow(kpc,-3)/4./Pi*1.e4;
00331   normel = norm;
00332 #else
00333   
00334   TParticle* protons = NULL;
00335   TParticle* electrons = NULL;
00336   if (prop_extracomp) electrons = particles.front();
00337   else {
00338     for (vector<TParticle*>::reverse_iterator ripart = particles.rbegin(); ripart != particles.rend(); ++ripart) {
00339       if ( (*ripart)->GetUid() == 1001 && !(*ripart)->GetIsSec() ) {
00340         protons = (*ripart);
00341       }
00342       if ((*ripart)->GetUid() == -1000 && !(*ripart)->GetIsSec()) electrons = (*ripart);
00343     }
00344   }
00345     
00346   if (protons) norm = protons->FindNormalization();
00347   if (electrons) normel = electrons->FindNormalization();
00348   
00349 #endif
00350 
00356   for (vector<TParticle*>::reverse_iterator ripart = particles.rbegin(); ripart != particles.rend(); ++ripart) { 
00357     if ( ripart != particles.rbegin() && ( (*ripart)->GetUid() == -1000 && !(*ripart)->GetIsSec() ) ) {
00358       if (output_ptr) (*ripart)->Print(output_ptr, normel);
00359       if (output_ptr_sp) (*ripart)->PrintSpectrum(output_ptr_sp, normel);
00360     }
00361     else {
00362       if (output_ptr) (*ripart)->Print(output_ptr, norm); 
00363       if (output_ptr_sp) (*ripart)->PrintSpectrum(output_ptr_sp, norm); 
00364     }
00365   }
00366 
00367   delete spallnet;
00368   delete xsecmodel;
00369   spallnet = NULL;
00370   xsecmodel = NULL;
00371   return ;
00372 }

Member Data Documentation

vector<TCREvolutorBasis*> DRAGON::alg [protected]

Container of possibly many resolution algorithms to be used in cascade.

Definition at line 59 of file dragon.h.

Referenced by DRAGON(), Run(), and ~DRAGON().

Galaxy* DRAGON::gal [protected]

Pointer to a Galaxy object, where all the ambient information is stored.

Definition at line 57 of file dragon.h.

Referenced by DRAGON(), particle_modulated_spectrum(), PrintChi2(), PrintRatios(), Run(), and ~DRAGON().

Input* DRAGON::in [protected]

Pointer to a Input object, an interface to user input.

Definition at line 58 of file dragon.h.

Referenced by DRAGON(), Print(), PrintChi2(), PrintRatios(), Run(), and ~DRAGON().

TNucleiList* DRAGON::list [protected]

List of nuclei to be propagated.

Definition at line 60 of file dragon.h.

Referenced by DRAGON(), Run(), and ~DRAGON().

fitsfile* DRAGON::output_ptr [protected]

Output structure for spectra and galactic densities of CRs. Used if fullstore == true.

See also:
fullstore

Definition at line 62 of file dragon.h.

Referenced by DRAGON(), Print(), Run(), and ~DRAGON().

fitsfile* DRAGON::output_ptr_sp [protected]

Output structure only for solar CR spectra. Used if partialstore == true.

See also:
partialstore

Definition at line 63 of file dragon.h.

Referenced by DRAGON(), Print(), Run(), and ~DRAGON().

vector<TParticle*> DRAGON::particles [protected]

Container of nuclei densities and spectra.

Definition at line 61 of file dragon.h.

Referenced by FindParticle(), Run(), and ~DRAGON().

int DRAGON::status [protected]

FITS error status.

Definition at line 64 of file dragon.h.

Referenced by DRAGON(), Print(), and ~DRAGON().

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:57 2010 for DRAGON by  doxygen 1.6.3