/*
  MAGIC Source Simulator. The macro takes a given spectral shape, 
  and  using MAGIC performance in the range 40GeV-16TeV from 
  Aleksic, J., et al., 2016, Astroparticle Physics, 72, 76 
  (please cite this paper if you want to use this macro in a publication)
  estimates what kind of signal can be observed by MAGIC.
  It computes significances of each spectral point according to
  Li & Ma 1983, ApJ, 272, 317, Eq. 17 definition.

  Terminal Output:
  For each estimated energy range one gets the number of excess events,
  signal-to-background ratio, significance, and information if a given
  bin satisfies the conditions for the detection. 

  Plot Output:
  Spectral points following the assumed spectrum are shown, their 
  error bars reflect the performance of the instrument for such a source.
  Significance for each bin is given (bins without detection
  have gray numbers. In the top of the plot a simple number using 
  information from all the shown bins is given to evaluate 
  if the source can be detected (roughly if sum of sigma_i / sqrt(N) >~5) 

  requirements:
  - ROOT (https://root.cern.ch), the macro was tested on 5.34/21 version

  caveats:
  - the macro is operating on estimated energy doing simple
  comparisons with differential in estimated energy rates seen for
  Crab Nebula, for soft sources the differences in energy migration 
  will result in different performance than calculated here
  - the threatment of extended sources is very approximate, only
  increase in background is taken into account (without energy
  dependence of PSF) for extended sources for extension >~0.4 deg the
  dependence on the offset from the centre of the camera will further
  worsen the performance w.r.t. the one calculated here
  - significances are given for each differential energy bin
  separatelly, but to detect a source one normally applies a cut that
  keeps a broad range of energies inside resulting in better integral
  sensitivity than differential one. Also, optimization of cuts for 
  a broad energy range usually results in somewhat better sensitivity
  than what one can get by simply integrating the used here signal in 
  differential energy bins. As a very crude approximation for detection
  capability we calculate here also a sum of significances of all the
  points in SED divided by the sqrt of the number of those points

  usage:
  - modify the "basic settings" below and run in terminal:
  root mss.cpp

  for questions ask Julian Sitarek (jsitarek [at] uni.lodz.pl)

  version history:
  1.1 (2017/10/16): fixed a bug which made the errorbars too small
 */

#include <TCanvas.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <TH1.h>
#include <TLatex.h>
#include <TLegend.h>
#include <TMath.h>

#include <iostream>

using namespace std;

// basic settings
const Double_t timeh=50; //[h], time of observations
const Double_t extension=0.; // [deg] extension radius of the source

// [cm^-2 s^-1 TeV^-1] assumed source spectrum (dN/dS dt dE) as a function of energy in GeV
// the string understands many mathematical expressions, some common examples below  
// const TString assumedstr="1.e-11 * pow(x/1000, -1.9)"; // simple power-law
// const TString assumedstr="1.e-11 * pow(x/1000, -1.9) * exp(-x/300)"; // power-law with an exponential cut-off
// const TString assumedstr="1.e-11 * pow(x/1000, -1.9) * exp(-pow(x/300, 0.5))"; // power-law with a non-exponential cut-off
// const TString assumedstr="3.39e-11 * pow(x/1000, -2.51-0.21*log10(x/1000))"; //Crab Nebula [Aleksic et al. 2016]
const TString assumedstr="0.1* exp(-x/200)* 3.39e-11 * pow(x/1000, -2.51-0.21*log10(x/1000))"; //10% Crab with 200 GeV cut-off

const TString srcname="SourceTest1"; // name of the source
const Bool_t ismidzd=kFALSE; // false: performance for 0-30deg zenith angles, true: for 30-45deg zenith
// you can check visibility of your source e.g. here: http://www.magic.iac.es/scheduler/

// advanced settings
const Int_t numoff=3; // number of background estimation regions
const Double_t minev=10; // minimum number of events
const Double_t minSBR=0.05; // minimum ratio of excess to background
const Double_t PSF=0.1; //[deg] PSF (for worsening the performance for extended sources)
const TString crabstr="3.39e-11 * pow(x/1000, -2.51-0.21*log10(x/1000))"; //Crab Nebula [Aleksic et al. 2016]
const Double_t eplotmin=40; //[GeV] x plot range
const Double_t eplotmax=20.e3; //[GeV] x plot range
const Double_t yplotmin=1.e-14; // [TeV cm^-2 s^-1] y plot range
const Double_t yplotmax=1.e-9; //[TeV cm^-2 s^-1] y plot range
const Double_t minerror=2; // showing only points with value > minerror * error
const Bool_t drawsigma=kTRUE; // whether to draw also sigmas on the plot

// global variables (DO NOT MODIFY)
const Int_t npoints = 13;
Double_t crabrate[npoints];// [min^-1]
Double_t bgdrate[npoints];// [min^-1]
Double_t enbins[npoints+1]; // [GeV]

void Prepare()
{
  for (int i=0; i<=npoints; i++)
    enbins[i]=100* pow(10, (i-2)*0.2);

  if (ismidzd)
    {
      crabrate[0]=0; // below threshold !!!
      crabrate[1]=0.404836; 
      crabrate[2]=3.17608; 
      crabrate[3]=2.67108; 
      crabrate[4]=2.86307; 
      crabrate[5]=1.76124;
      crabrate[6]=1.43988; 
      crabrate[7]=0.944385; 
      crabrate[8]=0.673335; 
      crabrate[9]=0.316263; 
      crabrate[10]=0.200331; 
      crabrate[11]=0.0991222; 
      crabrate[12]=0.0289831;

      bgdrate[0]=1.67777; // below threashold
      bgdrate[1]=2.91732; 
      bgdrate[2]=2.89228; 
      bgdrate[3]=0.542563; 
      bgdrate[4]=0.30467; 
      bgdrate[5]=0.0876449; 
      bgdrate[6]=0.0375621; 
      bgdrate[7]=0.0197085; 
      bgdrate[8]=0.0111295; 
      bgdrate[9]=0.00927459; 
      bgdrate[10]=0.00417356; 
      bgdrate[11]=0.00521696; 
      bgdrate[12]=0.000231865;      
    }
  else
    {
      // 0-30 deg values
      crabrate[0]=0.818446;
      crabrate[1]=3.01248;
      crabrate[2]=4.29046;
      crabrate[3]=3.3699;
      crabrate[4]=1.36207;
      crabrate[5]=1.21791;
      crabrate[6]=0.880268;
      crabrate[7]=0.579754;
      crabrate[8]=0.299179;
      crabrate[9]=0.166192;
      crabrate[10]=0.0931911;
      crabrate[11]=0.059986;
      crabrate[12]=0.017854;
      
      bgdrate[0]=3.66424;
      bgdrate[1]=4.05919;
      bgdrate[2]=2.41479;
      bgdrate[3]=0.543629;
      bgdrate[4]=0.0660764;
      bgdrate[5]=0.0270313;
      bgdrate[6]=0.0132653;
      bgdrate[7]=0.00592351;
      bgdrate[8]=0.00266975;
      bgdrate[9]=0.00200231;
      bgdrate[10]=0.00141831;
      bgdrate[11]=0.00458864;
      bgdrate[12]=0.0016686;
    }
}

// Li & Ma eq 17. significance, function abridged from MARS 
Double_t SignificanceLiMa(Double_t s, Double_t b, Double_t alpha)
{
    if( b==0 && s==0)
      return 0;

    b = b==0 ? FLT_MIN:b; // Guarantee that a value is output even in the case b or s == 0
    s = s==0 ? FLT_MIN:s; //   (which is not less allowed, possible, or meaningful than
                          //    doing it with b,s = 0.001)
    const Double_t sum = s+b;

    if (sum<0 || alpha<=0)
        return -1;
    const Double_t l = s*TMath::Log(s/sum*(alpha+1)/alpha);
    const Double_t m = b*TMath::Log(b/sum*(alpha+1)      );
    return l+m<0 ? -1 : TMath::Sqrt((l+m)*2);
}

Bool_t Checks()
{
  if (extension>1)
    {
      cout<<"Extension comparable to the size of the MAGIC camera cannot be simulated"<<endl;
      return kFALSE;
    }
  if((numoff<=0)||(numoff>7))
    {
      cout<<"Number of OFF estimation regions must be in the range 1-7"<<endl;
      return kFALSE;
    }
  if ((extension >0.5)&&(numoff>1))
    {
      cout<<"For large source extensions 1 OFF estimation region (numoff) should be used"<<endl;
      return kFALSE;
    }

  return kTRUE;
}
void mss()
{
  if (!Checks())
    {
      cout<<"exiting"<<endl;
      return;
    }

  Prepare();
  // preparing reference SED graph;
  TF1 *ffluxcrab = new TF1 ("ffluxcrab", crabstr, eplotmin, eplotmax);
  TF1 *fsedcrab = new TF1 ("fsedcrab", Form("1.e-6*x*x*%s", crabstr.Data()), eplotmin, eplotmax);
  fsedcrab->SetLineColor(kGray+2);

  TF1 *ffluxass = new TF1 ("ffluxass", assumedstr, eplotmin, eplotmax);
  TF1 *fsedass = new TF1 ("fsedass", Form("1.e-6*x*x*%s", assumedstr.Data()), eplotmin, eplotmax);
  fsedass->SetLineColor(kGreen);

  TGraphErrors *gsed = new TGraphErrors();


  TCanvas *can = new TCanvas ("can", "", 640, 480);
  can->SetLogx(); can->SetLogy();
  can->SetGrid();
  TH1F *hframe = can->DrawFrame(eplotmin, yplotmin, eplotmax, yplotmax);
  hframe->GetXaxis()->SetTitle("E [GeV]");
  hframe->GetYaxis()->SetTitle("E^{2} dN/dE [TeV cm^{-2} s^{-1}]");
  hframe->GetYaxis()->SetTitleOffset(1.3);
  fsedcrab->Draw("same");
  fsedass->Draw("same");

  Double_t sig_comb = 0.0;
  Int_t sig_comb_N = 0;

  for (int i=0; i<npoints; i++)
    {
      Double_t dummy[3];
      Double_t intcrab=ffluxcrab->Integral(enbins[i], enbins[i+1], dummy, 1.e-15);
      Double_t intass=ffluxass->Integral(enbins[i], enbins[i+1], dummy, 1.e-15);
      Double_t noff = bgdrate[i]*timeh*60; // number of off events
      noff*=(PSF*PSF+extension*extension)/(PSF*PSF); // larger integration cut due to extension

      Double_t dnoff=sqrt(noff/numoff); // error on the number of off events (computed from numoff regions)
      Double_t nexc = crabrate[i]*timeh*60 * intass/intcrab; 
      Double_t dexc = sqrt(nexc + noff + dnoff*dnoff); // error on the number of excess events 
      
      Double_t sigma = 0; 
      // for tiny excesses (1.e-7 events) the function below can have numerical problems, and either way sigma should be 0 then
      if (nexc>0.01) 
	sigma = SignificanceLiMa(nexc+noff, noff*numoff, 1./numoff);
      Bool_t detect=kFALSE; 
      if ((sigma>=5.0)&&(nexc/noff>minSBR)&&(nexc>minev))
	detect=kTRUE;
      cout<<Form("%.1f-%.1f GeV: exc. = %.1f ev., SBR=%.2f%%, sigma = %.1f", 
		 enbins[i], enbins[i+1], nexc, 100.*nexc/noff, sigma)<<
	  (detect?" DETECTION":"")<<endl;
      if (nexc>minerror*dexc)
	{
	  Double_t en = sqrt(enbins[i]*enbins[i+1]);
	  Double_t sed = fsedass->Eval(en);
	  Double_t dsed = sed * dexc/nexc;
	  gsed->SetPoint(gsed->GetN(), en, sed);
	  gsed->SetPointError(gsed->GetN()-1, 0, dsed);
	  sig_comb += sigma;
	  sig_comb_N++;
	  if (drawsigma)
	    {
	      TLatex *text = new TLatex (en, 2*sed, Form("%.1f#sigma", sigma));
	      text->SetTextColor(detect?kBlack:kGray);
	      text->SetTextSizePixels(14);
	      text->SetTextAngle(90);
	      text->Draw();
	    }
	}
    }
  if (sig_comb_N>0)
    sig_comb = sig_comb/sqrt(float(sig_comb_N));
  cout<<Form("Combined significance (using the %d data points shown in the SED) = %.1f", 
	     sig_comb_N, sig_comb) << endl;
  // TLatex *combsig = new TLatex(eplotmin*1.5, yplotmin*2, Form("#Sigma#sigma_{i}/#sqrt{N}=%.1f", sig_comb_N, sig_comb));
  if (sig_comb<4)
    cout<<"The source probably will not be detected by MAGIC in "<<timeh<<" hours"<<endl;
  else if (sig_comb<6)
    cout<<"The source probably might be detected by MAGIC in "<<timeh<<" hours"<<endl;
  else
    cout<<"The source probably will be detected by MAGIC in "<<timeh<<" hours"<<endl;

  gsed->SetMarkerStyle(20);
  
  TLegend *leg = new TLegend (0.5, 0.8, 0.99, 0.99);
  leg->AddEntry(fsedass, Form("%s (Assumed)", srcname.Data()), "L");
  leg->AddEntry(gsed, Form("expected MAGIC SED (T_{obs} = %.1f h)", timeh), "PE");
  leg->AddEntry(fsedcrab, "Crab (Aleksic et al 2016)", "L");

  gsed->Draw("PE");
  leg->Draw();
  if (drawsigma)
    {
      TLatex *text = new TLatex();
      text->DrawLatexNDC(0.2, 0.93,  Form("#Sigma#sigma_{i} / #sqrt{%i} = %.1f#sigma", sig_comb_N, sig_comb));
    }
}