#ifndef MARS_MCascade
#define MARS_MCascade

#ifndef ROOT_TH1
#include <TH1.h>
#endif

#ifndef MARS_MPairProduction
#include "MPairProduction.h"
#endif

class TList;
class TFile;
class TBranch;

class MPhoton;
class MElectron;

class MCascade : public TObject
{
private:
    Double_t fSrcZ;
    Double_t fSrcR;

    Double_t fB;
    Double_t fBubbleZ;

    Double_t fRuntime;

    Double_t fEHi;
    Double_t fELo;

    Double_t fSpectralIndex;
    Double_t fDisplayIndex;

    Int_t    fNumBins;
    Int_t    fNumMaxInvCompton;
    Double_t fRatioInvCompton;

    // -------------------------------

    MPairProduction fPair;     //!
    TBranch *fBranchGammas;    //!
    TBranch *fBranchElectrons; //!

    TH1D fHist; //!

    Bool_t fIsBatch; //!

    Double_t GetEnergy();
    TFile *OpenFile(TString fFilename);
    void CloseFile(TFile *fFile);
    void ProcessElectron(MElectron &e, TList &fListGammas);
    Bool_t ProcessGamma(MPhoton &p, Double_t weight, TList &fListElectrons);
    Double_t ProcessGammas(TList &fListGammas, TList &fListElectrons, Double_t weight);
    Double_t ProcessElectrons(TList &fListElectrons, TList &fListGammas);
    void ProcessPrimaryGamma(Double_t E, Double_t weight);

public:
    MCascade();
    ~MCascade();

    void SetSourceZ(Double_t z);
    void SetSourceRZ(Double_t r); // [kpc]
    void SetB(Double_t B) { fB = B*1e-4; } // [G=1e-4T] mean magnetic field
    void SetRuntime(Double_t min) { fRuntime = min*60; } // [s] maximum time to run the simulation
    void SetEnergyBins(Int_t n, Double_t lo, Double_t hi);
    void SetMaxInvCompton(Int_t max) { fNumMaxInvCompton = max; }  // maximum number of inv. Compton (-1 means infinite)
    void SetRatioInvCompton(Double_t r) { fRatioInvCompton=r; }
    void SetSpectralIndex(Double_t idx) { fSpectralIndex=idx; }
    void SetDisplayIndex(Double_t idx) { fDisplayIndex=idx; }
    void SetBradius(Double_t r); // [Mpc]

    void Run(TString filename, Bool_t draw);

    ClassDef(MCascade, 1)
};

#endif