// Ulli Schwanke, Louise Oakes
const double sv = 3e-26; //annihlation cross-section [cm^3/s]
const double kpc2cm = 3.08e21; //1 kpc = 3.08e21 cm
const double M = 500; //mass of the DM particle [GeV]
const double Br = 1e-4; //branching ratio

double rate(char* place,double radius,Double_t density,double sigma=sv){
  double V=4./3.*TMath::Pi()*pow(radius,3);
  double n = density/M;
  double r = sigma*V*n*n/2.0; //need to divide by 4.0 (instead of 2.0) if Dirac and not Majorana
  std::cout << place << " DM density " << density << " GeV/cm^3" << std::endl;
  std::cout << place << " annihilation rate is " << r << " Hz = " << r*3600*24 << " per day" << std::endl;
  return r;
}

void rates(){
  std::cout << "Annihilation cross section: " << sv << " cm^3/s" << std::endl;
  std::cout << "Mass of the DM particle   : " << M << " GeV" << std::endl;
  rate("Earth: ",6371e5,0.3); //radius 6371 km (in cm), 0.3 GeV/cm^3
  double r = rate("Toy  : ",0.2*kpc2cm,150); //radius 0.2kpc (in cm), 30 GeV/cm^3
  double flux = r*2.0*Br/4./TMath::Pi()/pow(8.5*kpc2cm,2);
  std::cout << "Toy  :  Flux is " << flux << " cm^-2 s^-1" << std::endl;
  double fluxCrab = 1.74e-10*0.2; //multiply 1.74e-10 cm-2 s-1 TeV-1 by 0.2 TeV = 200 GeV 
  std::cout << "Crab :  Flux is " << fluxCrab << " cm^-2 s^-1" << std::endl;
}